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I .   INTRODUCTION 

Ocean  Thermal  Energy  Conversion  (OTEC)  plants  will  be 
either  moored  or  dynamically  positioned  in  the  ocean  and, 
accordingly,  will  be  subject  to  wind,  wave  and  current  load- 
ing.  This  report  deals  with  the  motion  of  a  floating  body 
of  arbitrary  shape  subject  to  wave  motion.   The  interaction 
of  the  waves  with  the  floating  body  is  based  on  the  inviscid 
(potential  flow)  theory.   This  assumption  appears  to  be 
valid  as  long  as  the  body  involved  is  large  compared  to  the 
amplitude  of  the  relative  fluid  motion. 

This  report  describes  the  first  phase  of  a  project  to 
analytically  determine  the  dynamic  response  of  large  float- 
ing OTEC  plants  to  ocean  waves.   The  long  range  goal  of  the 
project  is  to  develop  analysis  and  computer  code  for  compu- 
tation of  the  dynamic  response  and  structural  loading  of 
OTEC  structures  of  rather  general  configuration.   Such 
configurations  are  considered  to  be  composed  of  a  large 
displacement  parts  where  diffraction  theory  will  be  required 
and/or  small  appendages  or  members  which  can  be  dealt  with 
by  use  of  a  Morison  equation  type  analysis.   This  phase  of 
the  work  deals  with  the  analysis  of  large  displacement 
bodies  of  arbitrary  shape.   The  analysis  is  based  on  the 
use  of  a  surface  distribution  of  sources  using  quadralateral 
source  patches  or  panels.   The  source  strengths  over  each 


panel  is  assumed  constant  in  this  phase  of  the  project.   The 
next  phase  of  the  work  will  investigate  the  use  of  triangular 
panels  where  the  source  strength  varies  linearly  between  the 
values  at  the  corner  node  points.   A  comparison  of  the  tri- 
angular patch  method  will  be  made  with  the  quadralateral 
source  patch  method  with  constant  source  strengths  over  the 
patch.   The  better  of  the  two  methods  will  be  used  in  the 
future  phases  of  this  work. 

In  particular,  this  report  deals  with  the  analytical 
evaluation  of  the  wave  forces  and  overturning  moment  acting 
on  a  large  bodyof  arbitrary  shapewithout  appendages  in  water 
of  finite  depth.   The  same  structure  is  also  considered  to 
oscillate  in  surge,  heave,  sway,  roll,  yaw  and  pitch  and  for 
such  motion  the  added  mass  and  damping  coefficients  are  deter' 
mined.   Then,  the  equations  of  motion  for  a  free-floating 
body  are  developed,  and  using  the  computed  wave  excitation 
forces  and  moments  and  added  mass  and  damping  coefficients 
the  response  of  the  floating  structure  is  evaluated. 

Although  the  problems  of  wave  interaction  with  a  fixed 
body  and  the  oscillation  of  the  same  body  in  otherwise 
still  fluid  appear  to  be  physically  distinct,  they  are 
mathematically  similar  and  as  a  result  are  dealt  with 
simultaneously  herein.   The  primary  mathematical  difference 
in  these  problems  is  the  kinematic  boundary  condition  which 
is  applied  on  the  immersed  surface.   In  all  cases,  however, 


the  velocity  of  the  fluid  in  the  normal  direction  relative 

to  the  immersed  surface  must  be  zero  but  this  statement  takes 

a  different  form  in  each  case. 

The  method  used  to  describe  the  fluid  motion  involves 
the  use  of  a  Green's  function  or  source  distribution.   The 
rigid,  immersed  surface  of  the  structure  is  represented  by  a 
surface  distribution  of  sources  which  is   assumed  to  be  con- 
stant over  quadrala teral  panels  and  the  zero  relative  normal 
velocity  condition  is  applied  in  order  to  determine  the 
strengths  of  the  sources.   Once  the  source  strength  distri- 
bution is,  known  the  velocity  potential  at  some  point  in  the 
fluid  region  is  determined  by  summing  the  effect  of  all  of 
the  sources  at  the  point  in  question. 

Several  authors  have  computed  hydrodynamic  coefficients 
associated  with  specific  shapes.   The  added  mass  and  damping 
coefficients  for  a  semi-immersed  floating  sphere  undergoing 
heaving  motion  has  been  obtained  by  Havelock  (2).   The  vertical 
wave  excitation  force  was  not  calculated  by  Havelock  but 
using  the  Haskind's  Relations  [Ref.3]  this  quantity  may  be 
determined  from  the  damping  coefficient.   Kim  (4)  determined 
the  added-mass  and  damping  coefficients  for  this  same  config- 
uration to  wave  excitation  in  heave  and  surge.   Although 
calculated  for  the  infinite  depth  case,  these  results  pro- 
vide an  important  comparison  for  establishing  confidence  in 
the  mathematical  validity  of  the  response  calculations  based 
on  the  present  analysis. 
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Several  authors  have  also  considered  the  wave  interaction 
with  large  displacement  fixed  and  floating  bodies.   In  the 
case  of  large  North  Sea  Gravity  platform  [Refs.  10  and  20] 
the  agreement  between  the  theory  and  experiment  is  almost  exact. 
In  such  structures  the  caisson  was  well-submerged  and  the 
superstructure  only  pierced  the  free-surface.   Thus,  nonlinear 
effects  which  tend  to  be  most  pronounced  near  the  free  surface 
were  absent.   Additional  experimental  data  and  comparisons 
with  the  linear  theory  for  the  case  of  fixed  bodies  which 
are  considered  to  be  significant  include  those  presented  by 
Garrison  and  Seetharama  Rao  (6)  for  the  case  of  a  bottom 
mounted  hemisphere.   Unfortunately,  the  wave  amplitudes  were 
rather  small.   Also,  Garrison  and  Chow  (7)  have  presented  a 
comparison  of  the  theoretical  results  based  on  the  distributed 
source  theory  similar  to  that  developed  herein  with  two  slightly 
different  submerged  oil  storage  tank  models.   Unfortunately, 
the  vertical  force  part  of  this  data  was  rather  poor  because 
the  vertical  model  supports  had  very  low  spring  rates  making 
the  natural  frequency  of  the  model  support  system  too  close  to 
the  wave  excitation  frequency.   It  is  significant,  however, 
that  these  tests  represented  North  Sea  design  conditions  and 
in  the  case  of  the  horizontal  force  quite  good  agreement  was 
obtained  between  the  linear  theory  and  experimental  results. 
A  third  experimental  study  involving  a  short  vertical  circular 
cylinder  (1.16  cylinder  radii  water  depth)  extending  from  the 
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bottom  through  the  free  surface  should  be  noted.   Chakrabarti  (8) 
has  made  measurements  of  the  horizontal  force  and  overturning 
moment  for  this  configuration  and  found  the  results  to  compare 
well  with  MacCamy  and  Fuch's  (1)  corresponding  theoretical  re- 
sults.  However,  Chakrabarti  gives  no  indication  of  the  wave 
height  involved  in  Ref .  (8) .   It  is  noted  that  the  MacCamy  and 
Fuchs  solution  is  simply  a  special  case  of  the  present  more 
general  method  so  that  agreement  with  this  solution  implies 
agreement  with  the  present  method. 

A  series  of  tests  have  been  conducted  by  Hogben  and 
Standing  (9)  of  the  National  Physical  Laboratory,  London, 
England  with  fixed  vertical  cylinders  of  square  and  circular 
cross-section.   They  also  made  calculations  using  a  computer 
program  based  on  diffraction  theory  similar  to  the  one  pre- 
sented here  and  good  agreement  between  theory  and  experiment 
has  been  observed. 

Only  limited  experimental  data  has  been  presented  in  the 
literature  for  the  dynamic  response  of  three-dimensional 
floating  structures  to  waves.   Among  these  the  results  of 
Faltinsen  and  Michelsen  (11)  are  the  most  complete.   They 
compared  calculations  of  diffraction  theory  with  the  heave 
response  of  a  floating  box  and  the  agreement  obtained  was 
quite  good.   Also  Garrison  (12)  showed  excellent  agreement 
between  the  results  of  diffraction  theory  and  the  heave  and 
pitch  response  of  a  disc  buoy.   However,  as  a  general  assess- 
ment of  the  experimental  data  available  in  the  literature, 
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for  the  hydrodynami c  coefficients  and  dynamic  response  of 
floating  bodies,  it  appears  that  the  amount  is  very  small  and 
the  understanding  of  viscous  effects  and  differences  or  sim- 
ilarities of  the  experimental  results  with  the  theory  is  non- 
existent.  This  is  an  area  where  additional  work  is  needed. 
2.   FORMULATION  OF  PROBLEM 

Consider  a  rigid  object  of  arbitrary  shape  and  having 
characteristic  lineal  dimension  a  with  center  of  gravity 
submerged  to  a  depth  d  beneath  the  free  surface  in  water  of 
depth  h  as  shown  in  Fig.  1.   The  structure  is  considered  to 
be  smooth  to  the  extent  that  its  unit  normal  vector  is  a 
continuous  function  and  it  may  or  may  not  intersect  the  bottom 
or  free  surface.   Two  coordinate  systems  are  identified,  x,y,z 

coordinates  with  origin  fixed  at  the  free  surface  and  the  body 

_  i  _  i  _  i  _     _ 

coordinates  x,y,z   coordinates  located  at  depth  y  =  -d.   The 

bars  over  the  symbols  denote  dimensional  quantities. 

The  mathematical  problem  which  is  now  established  is  that 
associated  with  the  fluid  motion,  pressures  and  resulting  forces 
induced  by  the  small  amplitude  oscillation  of  the  object  in 
its  six  degrees  of  freedom  as  well  as  the  fluid  motion  asso- 
ciated with  the  interaction  of  the  fixed  object  with  a  train 
of  regular  waves.   The  small  amplitude  oscillatory  motion  of 
the  structure  about  its  equilibrium  position  with  frequency  a 
is  described  by  the  relationship, 

X.(t)  =  X?  R  [e"i0t],  i  =  1,2,3  (2.1a) 

i         l   e 
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FIGURE  1   DEFINITIONS 
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O.(t)  =  0°  Re[e  iat],  i  =  4,5,6 


(2.1b) 


where  the  subscript  i  =  1,2,3  denotes  lineal  oscillation  in 

_  i  _  i  _  i 
the  x,y,z   directions,  respectively,  and  i  =  4,5,6  denote  angular 

oscillation  about  the  x,y,z   axes,  respectively.   The  real 


numbers   X?  and  J?  denote  the  amplitudes  of  the  motion.   The 

displacements,  X..  ,  X„  and  X   are  normally  referred  to  as  surge, 

heave  and  sway,  respectively,  while  the  angular  components  0, ,  0 

and  0,  are  called  roll,  yaw  and  pitch,  respectively.   The  second 
o 

problem  dealt  with  simultaneously  is  the  interaction  of  a  train 
of  regular  surface  waves  with  the  object  fixed  in  space.   The 
regular  incident  waves  of  wave  height  2^°  =  H  and  wave  length 
L  are  assumed  to  progress  in  the  positive  x-direction  and 
interact  with  the  fixed  structure. 

The  fluid  is  assumed  to  be  incompressible,  and  the  motion 
irrotational  and  harmonic  with  time  dependence   e      in  all 
cases.   It  follows,  therefore,  that  a  velocity  potential  exists 
such  that  the  fluid  velocity  vector  may  be  defined  as 


q   =  R  [V<J,  (x,y,z)e  10t],  j  =  1,2,  ..6 
J      "J 


(2.2) 


where   R  [  <t>  .  e     1  denotes  the  velocity  potential  associated 

e   J 

with  the  motion  induced  by  oscillations  in  the  six  degrees  of 

-    -f    9     ->-  3     ->  3  -*■-*■   f- 

freedom  and  V  =  i  +  j  —3  +  k  —3  ,  the  symbols  i,j,k  denoting 

3  x       3y       3  z 
the  unit  vectors  in  the  x,y,z  directions,  respectively. 

For  the  case  of  regular  wave  interaction  with  the  fixed 
structure,  the  velocity  potential  may  be  written  as  the  sum 

(j),=<l>o  +  <l>7  (23) 
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where  <$>  0  denotes  the  potential  associated  with  the  incident 
wave  in  the  absence  of  the  structure  and  cj>  7  denotes  the 
scattering  potential  which  results  from  the  presence  of  the 
structure  in  the  wave  train.   In  this  case,  the  fluid  velocity 
vector  is  given  by 


q»  =  Re[V(4>0  +  4»7)e  1Qt] 


(2.4) 


The  continuity  equation  shows  that  <J>  .  ,  (j  =  0,1,2,  ...7) 


must  satisfy  the  Laplace  equation 


^2 


<j>j  (x,y,z)  =  0 


(2.5) 


and  from  the  linearized  form  of  Bernoulli's  equation,  which 
is  applied  throughout,  the  dynamic  fluid  pressure  is  given  by, 


P   =  R  [i  pa  cf>   e  i0t],   j  =  1,2,  ..  .6 
J     ®        J 


(2.6) 


For  the  second  problem  involving  wave  interaction  with 
the  fixed  object,  the  pressure  is  given  by 


P'     =    Re[i    po(4)o    +    <f>7)     e    iat] 


(2.7) 


The  velocity  potentials  must  satisfy  certain  boundary  con- 
ditions in  addition  to  Eq.  (2.5).   These  include  the  linearized 
free  surface  boundary  condition, 

2 


3<l>  • 

— -J(x,0,z) 

3y 


g    J 


(x, 0,z)  =  0,  j  =  0,1,2,  ..  .7 


(2.8) 
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Equation  (2.8)  is  familiar  from  linear  wave  theory.  Also,  <p  . 

must  satisfy  the  kinematic  boundary  condition  on  the  bottom  in 
all  cases, 

li  (x,-h,i)  =0,   j  =  0,1,2,  ...7  (2.9) 

9y 

The  surface  of  the  structure  in  its  mean  position  is 
described  by 

S(x,y,i)  =  0  (2.10) 

When  the  structure  oscillates,  the  velocity  of  the  fluid  in 
the  direction  normal  to  the  immersed  surface  must  equal  the 
velocity  of  the  surface  normal  to  itself.   Also,  for  the  case 
of  wave  interaction  with  the  fixed  body  it  is  necessary  that 
the  normal  velocity  component  be  zero.   These  conditions  are 
stated  mathematically  as 


All     _ 
3h 

3<t>  2     _ 
35 

jH.3     = 
8n 

aj^   _ 

3n 

il5    - 

Ttn 

lie  = 

3n 

3i_7     _ 
3n 


-  ioXn  n 
1    x 


-ioX„n 
2    y 


iaX0n 
3    z 


-ioQ , [ (d+y ) n    - zn 
4  z         v 


-iaO  c [ zn       -    xn 
5  x  z 


-iaO,[xn       -     (d+y)n    ] 
6         y  x 

"  3n 


(a) 
(b) 
(c) 
(d) 
(e) 
(f) 
(g) 


(2.11) 
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where  n  =  in   +  in   +  kn   denotes  the  unit  normal  vector  on 
x      y       z 

the  surface  of  the  object  directed  outward  into  the  fluid  and 
d  denotes  the  depth  of  submergence  of  the  x,'y,'z 'coordinate 
origin.   Finally,  the  velocity  potentials  must  satisfy  the 
radiation  condition  which  allows  only  outgoing  waves, 

<K(r  9,y)-A  ±(®)\~li    cosh[k(y  +  h_Ueiki      -         _  (2.12) 

cosh (kh)  1 

-2-2^  -1  -  - 

where  r   =  [x   +  z  ]    and  6  =  tan   (z/x).   The  wave  number  is 

defined  as   k  =  2tt/L   where   L   denotes  the  wave  length  and 

is  related  to  the  frequency  of  the  motion  according  tc  the 


well-known  relationship, 

2 


-   =  k  tanh  (kh) 


(2.13) 


The  velocity  potential  of  the  incident  wave  alone  progres- 
sing in  the  positive  x  direction  which  satisfies  Eqs .  (2.5), 
(2.8)  and  (2.9)  is  given  by 


*o(x,y)     = 


•     n°  un   ruj.~\  l     ik(x    cos     ty   +    z    sin    ifO 

_    ±£_D*    cosh[k(h+y)  ]e  (2    14) 

a       cosh  (kE) 


where  as  indicated  in  Figure  1,   ^   denotes  the  angle  of  the 
incident  wave  measured  clockwise  from  the  x-axis  and   n^  =  H/2 
denotes  the  amplitude  of  the  incident  wave,  H  being  the  wave 
height.   Moreover,  from  the  linearized  form  of  Bernoulli's 
equation,  the  free  surface  of  the  incident  wave  alone  is 
evaluated  from 

(2.15) 


Using     (2.14)     this     gives, 


1A    =     i^    cos(kx    cos    ty    +    kz    sin    \p) 


(2.16) 


Thus,  for  purposes  of  reckoning  phase  angle  it  is  noted  that 
at  t  =  0  the  crest  of  the  incident  wave  is  just  passing  the 
coordinate  origin. 

For  convenience  in  carrying  out  the  solution  for  the  seven 
potentials,  <(> .  ,  and  to  show  clearly  the  dependence  of  the 
solution  on  the  parameter,   a  =  2TTa/L  =  ka,  the  relative  water 
depth,  h  =  h/a  and  the  relative  depth  of  submergence,  d  =  d/a, 
the  space  variables  and  amplitudes  are  first  made  dimens ionless 
with  the  characteristic  linear  dimension  of  the  object,  a, 


x  =  x/a,  y  =  y/a,  z  =  z/a,  r  =  r/a 

o      _  o  o        o 

Xi  =  X±/a,     (i  =  1,2,3),  Xi  =  9  ,  (i  =  4,5,6) 

°      — °  9  _ 

n*  =  n  /a  =  H/2a,  r   =  r../a,  V  =  o  a/g 


(2.17) 


and  then  the  d imens ionles s  potential  functions  u.  are  introduced 

3 

as  , 


io  <j> .  (x  ,  y  ,  z)  /  gaX  .     =    a    tanh(ah)     u.(x,y,z),       j     =    1,2,  ..6 


ia  <j>7  (x,y  ,  z)  /ga  r^    =    -a    u^x.y.z) 

o 

ia  <j>0  (x,y)  /ga  n^    =    -a    u0(x,y) 


(2.18a) 
(2.18b) 
(2.18c) 


The  complex  dynamic  pressure  amplitude  can  now  be  written 
by  use  of  the  linearized  form  of  Bernouli's   equation  (2.6)  and 
(2.7) ,  as 
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p.  =  a  tanh(ah)  u.(x,y,z),   j  =  1,2, ...6  (2.19a) 

cosh [a (h+y ) ] e  y  -a    u  (  x  ,  y  ,  z )      (2.19b) 

cosh (ah) 


where  the  complex  amplitudes  of  the  pressure,  p   are  defined  as 
P. 


p  gaX 


^"5"  "  R  [p  •  (x,y  ,  z)  e     ]  ,   j  =  1,2,  ...6 


e   3 


(2.20a) 


-iat 


P  gan*     e 


=  R  [p  ,(x,y,z)e  "UL] 


(2. 20b) 


The  boundary  value  problem  which  describes  the  fluid  motion 
arising  from  the  oscillation  of  the  rigid  object  in  its  six 
degrees  of  freedom  as  well  as  the  scattering  of  the  incident 
wave  may  now  be  written  concisely  in  terms  of  dimens ionless 
parameters.   The  potentials  u.(x,y,z),  i  =  1,2, ...7,  continuous 
in  the  fluid  region  is  sought  such  that: 


V  u  (x,y, z)  =  0 

3_u  .  (x  ,  o  ,  z) -a  tanh(ah)  u.(x,o,z)  =  0 
3y3  J 

9u  .  (x , -h , z)  =  0 

9_u  (x,y,z)  =  g  (x,y,z)  on  S(x,y,z)  =  0 
9nJ  3 

/„  Q   \  -i  /a\  -12   cosh  [  a  (h+y)  ]  e    1 

u.(r,e,y)-A.(e)r     L  ;   [     i  +o ,  r  - 

j  l     1      cosh (ah)  1 


(2.21a) 
(2.21b) 

(2.21c) 
(2.21d) 

(2.21e) 
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The  g.(x,y,z)  functions  represent  the  prescribed  functions  which 
depend  on  the  mode  of  oscillation  (j  =  1,2, ...6),   and  j  =  7 
corresponds  to  scattering.   These  functions  represent  the  dimen- 
sionless  form  of  (2.11)  and  are  given  by 


81  =  nx' 


,0  =  n 

2    y 


3     z 


g4  -  (d+y)nz-zny,  g5  =  zn^xn^  g^    =  xny-(d+y)n. 


(2.21f) 


-?  =  w — TT  tn   sinh  [a  (h+y)  ]+  i  cosh[a(h+y)  ] 

7    cosh(ah)    y  J  l  J 


(ncos  i>  +    n       sin  ^)  1  e 
x  z 


ia(x  cos  i>  +    -  sin  ^  ) 


Equations  (2.21)  define  seven  boundary-value  problems  which 
correspond  to  oscillation  of  the  immersed  surface  in  each  of 
the  six  degrees  of  freedom  and  to  scattering  of  the  incident 
wave  by  the  fixed  body.   The  problem  statement  for  each  of  the 
potentials  is  the  same;  the  only  difference  lies  in  the  form  of 
the  boundary  condition  to  be  applied  on  the  immersed  surface. 

3.   REPRESENTATION  OF  THE  POTENTIAL 

The  boundary  value  problem  for  oscillation  in  the  six  degrees 
of  freedom  and  scattering  of  the  incident  wave  is  specified  in 
(2.21).   The  solution,  (i.e.),  the  function  u.(x,y,z),  may  be 
represented  by  use  of  a  Green's  function  having  the  physical 
interpretation  of  a  point  wave  source  of  unit  strength.   These 
sources  are  distributed  over  the  surface  of  the  object  according 
to  the  source  strength  function, f,  so  that  the  potential  at 
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some  point  (x,y,z)  within  the  fluid  region  is  given  by  the 
surface  integral 

Uj(x,y,z)  =  -^  JJ       fj  (?,n,  C)   G(x,y,z;C,n,0  dS      (3.1) 

where  (?,n>0  represents  a  point  on  the  surface  of  the  struc- 
ture, G  denotes  the  Green's  function  (or  source  potential)  and 

-  _  2 
dS  =  dS/a   denotes  the  dimens ionless  surface  area  element.   The 

Green's  function  is  defined,  therefore,  as  the  function  which 

satisfies 


V  G(x,y,z;?,n, C)  =  6(x-£)  6(y-n)  6(z-c) 


(3.2) 


as  well  as  the  boundary  conditions  (2.21b,  c,  and  e).   Such 


a  function  is  given  by  Wehausen  and  Laitone  (13)  as 

G(x,y,z;S,n,C)  =  |  +  G* 

where     G*  =  — , 


(3.3a) 


-yh 


+  2P.V.  f™    (y+v)e    cosh[y(n+h)]  cosh[y(y+h)]  Jo(yr)dM   (3-3b) 
J0  y  sinh  (yh)  -v  cosh  (yh) 

2   2 
.  2fT(a  -v  )  cosh  [a  (n+h)  ]  cosh  [a(y+h)]    jn(ar) 

2      2 
a  h  -  v  h  +  v 


R  =  [(x-O2  +  (y-n)2  +  (z-02]h 
R'=  [(x-U2  +  (y  +  2h+n)2  +  (z-c)2]^ 


(3.3c) 
(3.3d) 


r  =  [(x-C)   +  (z-c)  V 

2- 
v  =  a  a/g  =  a  tanh(ah) 


(3.3e) 
(3.3f) 
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and  P.V.  indicates  principal  value  of  the  integral.   An  alter- 
nate series  form  of  the  Green's  function  is  also  given  by 
Wehausen  and  Laitone  (13)  as 

2   2 
G(x,y  ,  z;  £,  x\  ,  O     =       yv    "a    cosh  [  a  (h+y  )  ]  cosh  [  a  (h  +  n  )  ]  [  Y  (ar)-iJ  (ar)] 

a  h-v  h+v  °        ° 


(y. 


+  4 


(3.4a) 


v  h-v 


cos  [u  (y+h)  ]  cos  [y.  (n+h)  ]K  (p.  r) 
k  k         o   k 


where  J   and  Y   denote,  respectively,  Bessel  functions  of  the 
o       o  r        J  ' 

first  and  second  kind  of  order  zero  and  K   denotes  the  modified 

o 

Bessel  function  of  the  second  kind  of  order  zero.   The  quanti- 
ties y,  are  the  real  positive  roots  of  the  equation 


u ,   tan(ij,h)  +  v  =  0 


(3.4b) 


The  Green's  functions  specified  in  (3.3)  and  (3.4)  are 
equivalent.   However,  for  purposes  of  numerical  evaluation  it 
is  found  that  the  series  formulation  given  in  (3.4)  can  be 
evaluated  more  accurately  and  with  much  less  expenditure  of 
computer  time  than  the  integral  form  when  (ar)  is  not  too  small. 
On  the  other  hand,  for  small  and  zero  values  of  (ar)  the  integral 
form  given  in  (3.3)  is  necessary.   The  integral  form  of  the 
Green's  function  as  given  in  Eq .  (3.3)  may  be  integrated 
directly  or,  it  may  be  rearranged  as  in  Appendix  C  so  as  to  be 
somewhat  better  conditioned  for  numerical  evaluation.   In  the 
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computer  program  when  the  integral  form  of  the  Green's  function 
is  required,  the  form  given  in  Appendix  C  is  actually  used  for 
numerical  evaluation.   It  appears  that  it  is  somewhat  less  time 
consuming  to  evaluate  than  the  form  given  by  Eq.  (3.3).   The 
potentials  u.  are  represented  by  (3.1)  and  with  the  introduc- 
tion of  the  Green's  functions  (3.3)  and  (3.4) (or  equivalent  form 
given  in  Appendix  C)  the  only  unknown  in  (3.1)  is  the  source 
strength  function  f..   To  evaluate  this  it  is  necessary  to 
take  the  derivative  of  (3.1)  in  the  direction  normal  to  the 
immersed  surface  and  then  apply  the  boundary  condition  (2.21d). 
This  results  in  the  following  integral  equation  from  which  f  . 
is  to  be  determined: 


1  ff  c     ft  x  3G(x 


,y  ,  z;  £,n , C)dS  =  g  (x,y,z),  j  =  1,2.  ..7 


(3 


where  3G/  3n  is  obtained  using  the  form  VG . n  where  VG  is  deter- 
mined by  straightforward  differentiation  of  (3.3)  or  (3.4). 
The  normal  derivative  of  (3.3)  is 


3_G 
9n 


-  -3  [nx(x-U  +  ny(y-n)  +  n^z-?)  ] 

-  £,  3[n  (x-O  +  n  (y  +  2h  +  n)  +  n  (z-c)] 
R  J   x  y  z 

2  P.V.  f°°    ^(^+v)e"M^oshfu  (u+h)  1 

'  J°       u  sinh(yh)  -  v  cosh(yh) 

rcosh[y(y+h) ]     J, (ur) 


(3.6) 


lVMr;  1 

[n     (x-C)+n    (z-cj-n       sinh [ u (y+h) ] J     ( ur ) 

x  z  y  J 
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-1 


2   2 
2Tra(a  -v  )  cosh  [  a  (  n  +  h)  ] 

2       2 

a  h  -  v  h  +  v 


cosh [a(h+y)  ] J  (ar)  /  v 
1      (  n^(x-U+n^(z-Q  j 


-n   sinh  [a(h+y)]J  (ar) 


and  the  normal  derivative  of  the  alternate  rearranged  form  of 
Eq  .  (3.3)  is  given  in  Appendix  C.   The  normal  derivative  of 
the  series  form  given  in  (3.4)  is 


2   2 
9G  . .  2tt(v  -a  )a  cosh[a(h  +  n)]  f  .,,,,,  nr,  ,   \      ■  i     r        \  i 
r—  -      p 5 sinh[a(h+y)][f  (ar)-iJ  (ar)]n 


1 

-[nx(x-0+nz(z-0]  cosh  [a(h+y)]±  [y  1  (  ar )  -  i  J±  (  ar )] 


-£ 


^k(Mk  +V  )C°S  [Mk(T?+h)  ] 


y2h  +  v2h  -  v 


sin[uk(y+h) ]Kq (ykr) n 


K  (y  r) 
+  cos[p  (y+h)]  -£ [n  (x-O  +  n^z-O] 


(3.7) 


where  n  ,  n   and  n   denote  the  components  of  the  outward  unit 
x    y        z 

normal  vector  on  the  immersed  surface. 

Returning  to  the  integral  equation  (3.5),  it  is  evident 

that  the  8G/3n  occurring  therein  as  given  by  (3.6)  is  singular 

3 
like  1/R   as  the  point  (  £  ,  n  ,  O  approaches  the  point  (x,y,z). 

Thus,  special  care  must  be  taken  when  integrating  this  term 

of  (3.6)  over  the  singularity.   Moreover,  the  kernel  of  (3.1) 

is  singular  like  1/R  at  this  same  point  and  likewise  special 

attention  must  be  given  to  this  situation.   Taking  the  integral 
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over  the  complete  surface  of  the  object  as  the  sum  of  two  parts, 
one  part  being  a  small  circle  of  radius  r   about  the  point 
(x,y,z),  as  indicated  in  Figure  2,  and  the  second  part  being 
the  remainder,  S',  of  the  surface,  we  may  write  (3.1)  as 


x,y,z;£,n,0  dS 


1   ffcfr  s       dG      ( 

+  4T7/S  '  fj^'n'^)  ^n"  (x»y»z»  €  ,  n  ,  C )  ds  =  g,(x,y,z) 


(3.8) 


where  I  denotes  the  area  inside  the  small  circle  of  radius  r  . 

o 

Now,  using  the  expression  for  3G/3n  as  obtained  from  (3.3a) 
and  keeping  in  mind  that  f.(£,n>C)  is  a  well-behaved  function, 
we  may  take  the  limiting  case  as  I,  or  equivalently ,  r   tends 
to  zero  to  obtain 


lim  1  .  ,      s   /T3G*  ,_ 


(3.9) 


+ 


+ 


lim  1   ff     f     ,  .  1G  (x,y,z;C,n, 

2^0  4tt  JJs'   ij^^'Z)     9n 


O     dS  =  g  (x,y ,z) 


Letting  the  point  (x,y,z)  lie  along  the  normal  at  a  distance 
e  above  the  surface  as  indicated  in  Figure  2,  it  is  evident 
that  the  first  integral  in  (3.8)  can  be  expressed  as 
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o 

CL 

cr 
< 

z> 
CD 

z 
I/) 


C\J 
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lim      1       .     ,  n /T  9     ,1^.  -,  •  1       *     (  ,3     fro    27irdr 


r    +0 
o 


lim 

r    -H 
o 

e  +  0 


(3.10) 


=  r0^°   [4ir  f j (x»y» z)   2tt   ( 


o 


-1) ]     =    -    -    f,(x,y,z) 


provided  r   is  chosen  to  be  small  enough  that  the  surface  area 
c  o  ° 

inside  r   may  be  considered  to  be  plane.   Moreover,  since  G* 
o 

occurring  in  (3.9)  is  regular,  the  second  integral  vanishes 
in  the  limit  as  Z+0  giving  the  integral  equation 

-f  (x,y,z)  +77j^fj(C,n,C)  ff  U'y'Z;  e.n.OdS  =  2  g.(x,y,z) 


(3.11) 


In  (3.11)  the  singular  point  at  (x,y,z)  is  considered  to  be 
excluded  from  the  surface  S. 

As  noted,  in  (3.1)  the  Green's  function  is  also  singular 
and  special  consideration  is  needed.   However,  using  a  procedure 
similar  to  the  above,  it  can  be  shown  that  this  1/R  singularity 
contributes  nothing  to  the  surface  integral  in  (3.1).   Thus, 
(3.1)  may  stand  as  is  without  modification  and  the  point  (x,y,z) 
considered  to  be  excluded  from  the  surface  S. 

However,  the  integral  of  1/R  over  the  facet  of  area  AS 
outside  the  singular  point  is  not  zero.   This  will  be  discussed 
in  detail  in  Appendix  B. 
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4.   HYDRODYNAMIC  FORCES  AND  MOMENTS 

The  forces  and  moments  caused  by  the  dynamic  fluid  pressure 
acting  upon  the  immersed  surface  of  the  structure  may  be  ob- 
tained from  the  integrals, 


F±j(t)  =  -(1.0  or  a)JJs?.gi    dS,  i,j  =  1,2,. 


(4.1a) 


F.(t)   =  -(1.0  or  iO^P'g.dS,  i  =  1,2, 


(4.1b) 


where  F. .(t)  denotes  the  i-th  component  of  load  arising  from 
the  j-th  component  of  motion  and  F.(t)  denotes  the  i-th  compon- 
ent of  wave  force  (or  moment).   The  coefficient  1.0  is  to  be 
used  in  the  case  of  a  force  (i  =  1,2,3)  while  a  is  to  be  used 
when  F  denotes  a  moment  (i  =  4,5,6) .   The  sign  convention  of 
the  forces,  moments,  displacements,  velocities,  etc.  follow 
the  right  hand  rule. 

It  is  convenient  and  conventional  to  place  the  forces  and 
moments  in  d imens ionless  coefficient  form.   Accordingly,  the 
following  complex  wave  force  (and  moment)  coefficients  are 
defined  using  (4.1b)  : 


F.,    ,el6i 
l (max) 

-3  0 

pga  n 


F.,    .ei6i 
l (max) 

-4  0 
pga  n 


,   i  =  1,2,3 


,   i  =  4,5,6 


(force) 


(moment) 


(4.2a) 


(4.2b) 
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The  symbol  F.(max)  denotes  the  maximum  value  of  the  oscillatory 

force  (or  moment)  and  is  taken  as  positive.   The  wave  force 

may,  therefore,  be  defined  by  the  complex  number  C.  or  by  the 

magnitude  of  C.  and  the  phase  angle  6..   Once  Ic.l  and  6.  are 
1  l  '  l  '        l 

known,  the  wave  force  (or  moment)  is  expressed  as  a  function 
of  time  as 

F.(t)  =  (1.0  or  a)p ga3(H/2a) | C. |  cos  (6.  -  ot)        (4.3) 

where  the  coefficient  (1.0)  is  applicable  in  the  case  of  a 
force  (i  =  1,2,3)  and  (a)  is  applied  when  (i  =  4,5,6)  and  F.(t) 
denotes  a  moment.   The  dimens ionless  wave  amplitude  is  defined 

o 

as  n^  =  H/2a.   It  is  further  recalled  that  all  phase  angles 
are  measured  in  relation  to  the  incident  wave;  at  t  =  0  the 
crest  of  the  undisturbed  incident  wave  is  located  at  the  coor- 
dinate origin . 

According  to  the  definitions  (4.2),  C.  may  be  expressed, 
using  (4.1b)  with  (2.20b)  and  (2.19a),  as 

ff  r  ,  .  cos  h  [  a  (h+y  )  ]       ia(x    cosifj+z    sinri/)  , 

C.  =  /  /_[a  u_(x,y,z)  -  £ — )  u\  e 

l    JJS  7  cosh  (ah) 


(4.4) 


.gjL(x,y,z)  dS,   i  =  1,2,  ...6 


Thus,  once  u7(x,y,z)  at  all  points  on  the  immersed  surface  is 

determined,  the  complex  wave  force  (or  moment)  coefficient  C. 

may  be  evaluated  by  evaluating  the  surface  integral  indicated 
in  (4.4) . 
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Turning  now  to  the  force  (or  moment)  produced  by  oscilla- 
tions in  the  various  degrees  of  freedom,  (2.19a)  and  (2.20a) 
may  be  used  in  conjunction  with  (4.1a)  to  obtain 

F±  (t)  — (1.0  or  a)pga3X°  vRg  [  e"10  tJJs   Uj  (  x  ,  y  ,  z)  g.  (  x  ,  y  ,  z  )  dS  ]  , 

i,j  =  1,2.  .6 


(4.5) 


where  as  previously  noted  (1.0)  corresponds  to  the  case  of  a 

force  (i  =  1,2,3)  and  (a)  corresponds  to  the  moment  (i  =  4,5,6)  . 

Furthermore,  if  the  two  dimens ionless  real  numbers  M. .  and  N. . 

ij        ij 


are  defined  as 


M.  .  =  -Re 


N 


jjs   uj (x,y , z) gi(x,y , z)  dS 
=  -Im /ys  u  (x,y  ,  z)  g±  (x,y  ,  z)  dS 


(4.6) 
(4.7) 


where  Re  and  Im  denote  real  part  and  imaginary  part,  respectively, 
then  (4.5)  may  be  written  as, 


F   (t)=(1.0  or  a)  p  ga3X%Re  [  (M.  .+±K±     )  e  iat],  i,j=l,2, 


(4.8) 


where,  as  in  (4.5),  the  coefficient  (1.0)  corresponds  to  i=l,2,3 
and  (a)  corresponds  to  1=4,5,6. 

Equation  (4.8)  may  be  further  rearranged  by  defining  the 
dimens ionless  parameters, 


and 


C  ,(t) 

iJ 


c   (t) 


F. . (t) 

_^J 

^3 


F.  ,(t) 

_4 
P  ga 


i  =  1,2,3;  j  =  1,2, 


4,5,6;  j  =  1,2,  ..  .6 


(4.9a) 


(4.9b) 
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2- 
With  (4.9)  and  using  the  definition  v  =o  a/g  ,(4.8)  may  be 

wr i  t ten , 


C . . ( t )  =  --M ..X.(t)-—  N..X.(t),  i,j  =  1,2, ...6 
ij         g  ij  3  g   ij  3 


(4.10) 


M..  and  N..  (i,j  =  1,2,  ...6)  represent  the  inertia  and 
damping  tensors,  respectively,  having  36  elements  each.   The 
first  index  denotes  the  direction  of  the  force  or  moment  and 
the  second  index  is  associated  with  the  component  of  the  motion 
It  can  be  shown,  moreover,  that  these  tensors  are  symmetrical 
so  that 


and 


M. .  =  M. . 


N  .  .  =  N  .  . 


(4.11a) 


(4.11b) 


Except  for  the  diagonal  terms  it  is  rather  difficult  to 
give  simple  definitions  for  the  dimens ionless  added  mass  and 
damping  tensors,  M. .  and  N. .,  respectively.   They  are  best 
defined  by  expressing  the  force  (or  moment)  acting  on  the  immersed 
surface  associated  with  the  dynamic  response  in  terms  of  these 
parameters.   Equation  (4.10)  gives: 


■«--{Y> 


r  4        4 

pa   M  .  .  X  .  +  pa  a   N 
i]  J. 


ij^l-  U-Aisie)  '  J=1'2"--6 


(4.12) 


where  when  (i  =  1,2,3)  F.  .  represents  a  force  and  when  (i  =  4,5,6) 

F.  .  represents  a  moment.   The  parameters  X.  and  X.  denote  the 
ij  3  3 

first  and  second  time  derivatives  of  the  dimens ionless  dis- 
placements : 
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X.  =  X./a      and       X.  =  X./a,    j  =  1,2,3 


(4.13a) 


X.  =  0  . 
3  3 


and 


Xj  =  0  ,       j  =  4,5,6 


(4.13b) 


5.   NUMERICAL  SOLUTION 

The  primary  objective  at  this  juncture  is  to  solve  the 
integral  equation  (3.11)  which  has  been  established  for  the 
source  strength  function,  f  .  (  £  ,  n  ,  L, )  .  Once  this  function  is 
obtained,  u.  can  be  determined  from  (3.1)  by  evaluating  the 
surface  integral.  Then,  all  other  physical  quantities  such 
as  pressure  and  resulting  forces  and  moments  may  be  determined. 

A  numerical  procedure  can  be  devised  by  approximating  the 
actual  immersed  surface  of  the  structure  by  a  contour  composed 
of  a  large  but  finite  number  of  facets.   In  the  limit  as  the 
number  of  facets  increases  and  the  size  decreases,  the  approx- 
imate contour  approaches  the  actual  contour.   Thus,  we  may  assume 
that  in  the  limit  as  the  number  of  facets  approaches  infinity, 
the  numerical  scheme  converges. 

To  begin,  the  immersed  surface  is  subdivided  as  indicated 
in  Figure  3  and  an  index  is  assigned  to  the  nodal  point  (centroid 
point)  of  each  of  the  small  facets.   It  is  required  that  (3.11) 
be  satisfied  at  each  of  these  nodal  points.   This  is  a  relax- 
ation of  the  requirement  implicit  in  (3.11)  that  the  equation 
be  satisfied  at  every  point  on  the  immersed  surface. 
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panel  node 
corner  node 


FIGURE  3   NUMERICAL  SCHEME 


3  4 


Recognizing  that  the  source  strength,  f.(£,n,C)>  is  a  regular 
well-behaved  function,  we  may  approximate  (3.11)  by 

1    N  ff         9G 

-fn(xi,yi,z.)  +  —  .E  fn(x.,y.,z.)yyAS_  —  (x-.y-.z.;  x,y,z)  dS 


2gn(x.,y.,zi) 


(5.1) 


where  AS.  denotes  the  area  of  the  i-th  facet  and  N  denotes  the 
J 

total  number  of  facets.   In  shorthand  notation,  (5.1)  may  be  written 


-f     +   f     a  .  . 
n .        n  .     ij 
J 


2g 


(5.2) 


where  n  =  1,  2,  ...7  corresponds  to  the  six  degrees  of  freedom  and 
scattering,  the  repeated  index  indicates  summation  as  usual  and  a.. 


is  def  ined  as 


a  .  . 


=  -h/Ss.       |f  (x-.Y-.z.;  x,y,z)  dS    (5.3) 


Eq.  (5.2)  is  a  complex  matrix  equation  which  may  be  solved  using  a 

digital  computer  to  obtain  the  source  strength  f    once  the  elements 

i 
of  the  square  matrix  a. .  are  calculated. 

As  an  approximation,  the  matrix  a. .  may  be  evaluated  by  using 

the  value  of  3G/3n  at  the  nodal  point  to  represent  the  mean-value 

over  the  elemental  surface  area  of  the  facet  and  thereby  further 

approximate  (5.3)  by 


ij 


1 
2tt 


—  (x±,y±,z±;     x,y,z)  AS. 


(5.4) 


where  3G/3n  may  be  evaluated  by  use  of  either  (3.6)  or  (3.7)  .   The 
choice  between  the  two  forms  depends  on  the  value  of  r.   When  r  is 
small,  the  infinite  integral  occurring  in  (3.6)  is  fairly  rapidly 
convergent  while  many  terms  of  the  series  indicated  in  (3.7)  are 
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required.   However,  when  r  is  large,  the  series  occurring  in  (3.7) 

converged  rather  rapidly  and,  consequently,  only  a  few  terms  are 

needed . 

The  potential  function  u.  may  be  obtained  from  (3.1)  once  the 

3 

source  strength  f .  is  determined  at  the  nodal  points.   In  a  manner 

similar  to  that  used  to  solve  the  integral  equation  (3.11),  the 

source  strength  is  taken  outside  the  surface  integral  and  (3.1)  is 
wr  it  t  en  as 

u^x^y^z.)  =  ±  fn(xj'yj'ZJ)i&S.G(xi'yi'ai;  X'Y'Z)  dS 

J=1  J  (5.5) 

where  as  in  (5.1)  n  =  1,2,  ...7,  and  the  integer  N  denotes  the  total 

number  of  nodal  points  on  the  immersed  surface.   In  indicial  notation 

(5.5)  becomes 

u     =   f     3  .  .  (5.6) 

n  .       n.    ij 

1        J 
where  the  square  matrix  3..  is  defined  as  the  integral  over  the 

panel  of  area  AS.  as 

J 

3ij       =      h      fLs       G^±>y±>*±->     x>y>z)     dS  (5.7) 

J 

Taking  the  value  of  G  at  the  nodal  point  to  represent  the  mean-value 
over  the  facet,  (5.7)  may  be  approximated  for  large  R  as 

B.j   =  1~       G(x.,y.,z.;  x.,y.,z.)  AS.  (5.8) 

where  G  is  evaluated  by  use  of  either  (3.3)  or  (3.4),  depending  on 
the  value  of  r.   When  R  is  not  large  (5.7)  must  be  used  and  the 
integration  of  1/R  term  in  G  is  discussed  in  Appendix  B. 

The  first  step  in  the  numerical  solution  to  the  problem  is  to 
evaluate  the  matrices,  a.,  and  3...   For  example,  the  point  (x,y,z) 
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on  the  immersed  surface  is  denoted  by  i  and  the  point  where  the 
source  is  located  (£,n,C)  is  denoted  by  j.   Thus,  using  the  series 
form  of  the  Green's  function  given  in  (3.4)  or  its  derivative  given 
in  (3.7),  there  is  no  difficulty  in  numerical  evaluation  once  the 
location  of  the  point  i  and  j  are  specified.   However,  the  integral 
form  of  G  and  its  derivatives  as  specified  by  (3.3)  and  (3.6)  re- 
quires some  special  consideration.   The  infinite  integral  which 
occurs  in  both  of  these  expressions  poses  two  difficulties  with 
respect  to  numerical  evaluation;  the  integral  has  an  infinite 

upper  limit  and  the  integrand  is  singular  at  y=y   where  y   is  the 

oo 


root  of 


u    tanh( y  h)  -  v 
o  o 


0 


(5.9) 


The  root  is  simply  equal  to  "a"  in  view  of  (3.3f)  .   However,  these 
difficulties  can  be  overcome  by  recognizing  that  the  integrand  is 
singular  like  l/(y-y  ),  subtracting  the  singularity  out  of  the 
integral  and  carrying  out  its  integration  analytically.   The  re- 
mainder of  the  integral  is  then  carried  out  numerically  and  the 
upper  limit  is  replaced  by  a  suitable  large  number  such  that  con- 
vergence is  assured. 

For  purposes  of  illustration  let  the  infinite  integral  in 
(3.3)  be  denoted  by  I  so  that 


I  =  P 


J   0   ( u-u  ) 


(5.10) 


wher  e 


Q(y)  = 


P(y)  (m-m   ) 

0 

y  tanh  ( yh)  -  v 


(5.11) 
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and  in  the  case  of  (3.3)  the  function  P(y)  is  defined  as 

n/  N     (y+v)e  J  cosh[p(h+n)]cosh[u(h+y)]  J  (yr)      ,             . 

P(y)  =  —, — r-r o          (J.I^J 

cosh (Mh) 


The  integral  in  (5.10)  may  be  rearranged  and  written  in  the 
f  o  rm 


.   /1Zyo     Q(y)-Q(yo)       dy    +    Q(y    )p>v  f2"o         ± 

J  °         (y-y J  °         J  °       (y-  v 


dy 


(5.13) 


/oo 
Q(y)  Hl, 

2y0    (y^lT)    dy 

where  the  principal  value  of  the  second  integral  can  easily 
be  shown  to  be  zero  and  the  other  two  integrals  are  proper 
and  easily  evaluated  by  numerical  integration.   Thus,  I  may 
be  evaluated  by  numerical  integration  of 


-/"  y(q(y)-Q(yn))  dy+  f       q(y)dy 


(5.14) 


where  Q  (y  )  represents  the  limit  of  (5.11)  as  ii^v       where 

y   =  a.   The  result  of  this  limiting  process  results  in  the 

following  expression: 

0(U     )     = P^n) 

VVMoy  tanh(y    h)|"l-y    h    tanh     (y    h)l+u    h  (5.15) 

O     L     O  O    -1     O 


or,  using  (3.3f)  and 


the  fact  that  y  =a,  (5.15)  becomes 


c\(       \  a  P(a) 

Q(yo}  =  v+(a2V)h 


(5.16) 


For  the  case  of  deep  water  (h->°°  )  we  find  v->a  so  that 

Q(yQ)  =  P(a) 


(5.17) 
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The  same  method  of  integration  may  be  used  to  evaluate  the 
infinite  integral  in  the  derivatives  of  G  as  well  as  in  G  itself 
The  only  difference  is  that  P ( u )  is  defined  differently  for  the 
other  cases  but  can  easily  be  determined  by  a  comparison  of 
(5.10)  and  (5.11)  with  the  particular  infinite  integral  to  be 
evalua  ted . 

An  alternate  method  of  evaluating  the  integral  form  of 
the  Green's  function  may  also  be  used.   This  procedure  as 
described  in  Appendix  C  involves  converting  the  infinite  upper 
limits  and  appears  to  be  particularly  useful  when  the  wave 
period  is  large.   The  computer  program  actually  uses  this  form 
for  numerical  calculation. 

There  is  one  further   significant  difficulty  in  evaluating 
a.  .,  and  3.  .•   When  point  i  is  distant  from  the  source  located 
at  point  j  it  is  adequate  to  use  approximates  such  as  Eq .  (5.4), 
and  (5.8) .   However,  when  i  =  j  or  point  i  is  near  point  j  the 
singular  term  in  the  Green's  function  which  is  of  the  form  1/R 
does  not  vary  slowly  over  the  panel  of  area  AS..   Thus,  when 
the  node  point  designated  by  the  index  i  is  rather  close  to 
the  j-th  node  point  the  1/R  term  in  G  and  derivatives  of  1/R 
which  occurs  in  a.  .,  must  be  integrated  over  AS.  rather  than 
approximated.   This  integration  is  described  in  detail  in 
Appendix  B. 
Specification  of  the  Subdivision  Scheme 

The  numerical  scheme  starts  with  specifying  panel  corner 
points  on  the  immersed  surface.   Each  corner  node  point  is 
assigned  an  index  and  then  a  correspondence  table  is  introduced 
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which  specifies  which  four  corner  node  points  make-up  a  given 
panel.   Thus,  an  index  is  assigned  to  each  panel. 

Given  the  coordinates  of  the  four  corners  of  a  given  panel 
it  is  rather  easy  to  compute  its  area  by  breaking  it  up  into 
two  triangular  areas.   Using  the  areas  and  centroids  of  the 
two  triangular  portions  of  the  polygon  the  centroid  of  the 
polygon  is  computed.   Finally,  the  unit  normal  vector  for  the 
panel  is  determined  by  taking  the  cross-product  of  vectors 
running  across  the  two  diagonals  of  the  polygon  panel.   These 
calculations  are  described  in  detail  in  Appendix  A. 

6  .   DYNAMIC  RESPONSE  FOR  A  FLOATING  BODY 

In  this  section  the  equations  of  motion  for  a  floating 
body  in  waves  are  developed.   These  equations  are  then  applied 
using  hydrodynamic  coefficients,  which  include  the  wave  excita- 
tion forces  and  moments  as  well  as  added  mass  and  damping  co- 
efficients, in  order  to  compute  the  dynamic  response. 

The  problem  under  consideration  is  represented  schematically 
in  Figure  1.   In  the  following,  the  equations  of  motion  are 
written  with  respect  to  the  center  of  gravity.   Thus,  the  origin 
of  the  body  coordinates  is  assumed  to  be  located  at  the  center 
of  gravity  of  the  floating  structure  and  d  is  defined  as  the 
d imens i onles s  depth  to  that  point. 

The  small  amplitude  displacement  of  the  body  center  of 
mass  with  respect  to  its  mean  position  in  the  inertial  reference 
frame  is  described  by  the  three  coordinates  X  (t),  X„(t)  and 
X^(t)  which  are  referred  to  as  surge,  heave,  and  sway,  respec- 
tively.  The  small  angular  displacements  of  the  body  about  the 
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x',  y'  and  z'   axes  are  denoted  by  0.,  0  ,  0   and  are  referred 

H  d  6 

to  as  roll,  yaw  and  pitch,  respectively. 

The  equations  of  motion  linearized  with  respect  to  the 
small  angular  displacements  of  the  body  may  now  be  written  as 
f o 1 lows : 


F1T(t)     =    mX1(t) 


F2     (t)     =    mX2(t) 


F3     (t)     =    mX3(t) 


F.      (t)     =     I     ,      ,     0.     -     I     ,      ,     0_        -        I    t »  q 
4  x    x         4  x    y  5  xz6 


(6.1a) 
(6.1b) 
(6.1c) 
(6. Id) 


F       (t)     =    I     ,     ,     0,    -    I  i  «    0,     -    I  ,   ,  0. 
5  y    y  5  y z       6  yx       4 


F ,     (t)     =    I     ,      ,0-Iti0.-I,      ,     0 
6  zz  6  zx4  zy  5 


(6  .le 
(6. If) 


T         T  T 

where  F..  (t),  F  „  (t)  and  F„  (t)  denote  the  three  components  of 

T         T 
the  total  external  force  acting  on  the  body  and  F   (t),  F    (t) 

T 
and  F,  (t)  denote  the  three  components  of  the  total  external 
6 

moment.   The  symbol  m  denotes  the  body  mass  which  equals  the 

displaced  mass.   The  moments  of  inertia  are  defined,  typically, 

as 

I  ,  ,  =  I  ,  ,  =  /-  x'y'  dm  (6.2) 

x  y      y  z     Jm  J 

where  the  integration  is  to  be  carried  out  over  the  complete 
mass  of  the  body.   For  bodies  having  symmetry  with  respect  to 
the  (x'-y1)  and  (y'-z1)  planes,  all  of  the  products  of  inertia 


41 


vanish.   Although  this  type  of  symmetry  is  common  to  most  ocean 
structures,  there  is  no  need  to  apply  the  limitation  at  this 
point  since  the  inclusion  of  the  product  of  inertia  terms  do 
not,  in  principle,  complicate  the  development. 

At  this  juncture  it  is  convenient  to  place  the  equations 
of  motion  in  dimens ionless  form  and  replace  the  x,  y,  z  sub- 
script system  with  a  less  cumbersome  indicial  system.   We  may 


define  the  following  dimens ionless  parameters  for  this  purpose 

_3 

m^_=Tn     =m     ssm     =tt. 

2  3     3  2 
.  -5  .  _5  _5 

m 


mll  ™  m22  =  m33  =  m/Pa  »  mi2  =  m21  =  m31  =  m13  = 


-5  ,  _5  _5 

l44  "  Xx'x'  /pa  »  m55  =  Iv'v'  /pa  '  m66  =  1i,i'  /pa 


m45=m54 


15 


--I-,-,  /pa  ,  m46=m6r^ti,  /pa  , 
X1(t)=X1(t)/I,  X2(t)=X2(t)/i,  X3(t)-X3(t)/a 

6 


-5  ,  _5 

m65=m56  =  -1y'it  /pa 


X4(t)=04(t),  X5(t)=05(t),  X6(t)=0£(t) 


f1T(t)=F1T(t)/pgi,  f2T(t)=F2T(t)/pga3,  f3T(t)=F3T(t)/pgi 
f4T(t)=F4T(t)/pSi  •  fcT(t)  =  F,T(t)/pgi4,  f,T(t)  =  F  T(t)/pgi4 


where,  as  previously  given,  a  denotes  the  characteristic  dimen- 
sion of  the  body,  p  denotes  the  fluid  density  and  g  denotes 
the  gravitational  constant. 
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Using  these  definitions,  (6.1)  may  be  condensed  to  the 


form 


f  .T(t)  =  -  m.  .  X.  (t) 
i         g   ij   J 


(6.4) 


where  ij  takes  on  values  1,2,  ...6,  and  the  repeated  index  denotes 
summation  as  usual. 

Equation  (6.4)  represents  the  six  equations  of  motion  with 

T 
f.  (t)  denoting  the  d imens ionles s  total  external  force  or  moment 

coefficients  as  defined  by  (6.3).   For  free-floating  bodies  these 

coefficients  represent  the  contributions  from  the  surrounding 

fluid  only  and  are  generally  considered  to  be  composed  of  three 

parts:     (a)  the  wave  excitation  forces  and  moments,  (b)  the 

dynamic  forces  and  moments  caused  by  the  motion  of  the  body, 

and  (c)  the  hydrostatic  forces  and  moments  caused  by  the  dis- 

placement  of  the  body.   For  moored  bodies  the  forces  and  moments 

associated  with  the  mooring  lines  must  be  included  and  for  the 

linear  problem  these  contributions  may  be  determined  separately 

and  superimposed. 

In  accordance  with  this  idea  the  three  contributions  to 

T 
the  force  (or  moment)  coefficient  f .  (t)  may  be  expressed  as 

the  linear  combination: 


,T(t)  =  c.(t)  +  Y      [c..(t)  +  k..(t)  +  k'    (t)] 

j-l 


(6.5) 
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wh  ere 


c.(t)  =  force  or  moment  coefficient  associated  with  wave 

excitation,  F.(t)/Pga   (i=l,2,3);  F.(t)/Pgl  (i=4,5,6) 
c. .(t)  =  force  or  moment  coefficient  associated  with  the 

dynamic  response  of  the  body,  see  Eqs.  (4.8-4.10). 
k. .(t)  =  force  or  moment  coefficient  associated  with  linear 
or  angular  displacement  which  results  from  hydro- 
static pressures, 
k'  .  .(t)  =  force  or  moment  coefficient  associated  with  linear 
or  angular  displacement  resulting  from  elastic 
constraints  (mooring  lines). 
The  force  ormoment    coefficient   c.(t)  may  be  expressed 
as 

c  (t)  =  n°  R  UcJe1  *  e"i0t],   i  =  1,2,  ..6     (6.6) 
J-         £   e     J. 

where  i  =  1,2,3  refers  to  the  three  components  of  the  force 
and  i  =  4,5,6  refers  to  the  three  components  of  the  moment. 

The  dimensionless  coefficients  c.(t)  are  defined  according  to 

T 
the  definition  of  f .  (t)  as  the  force  made  dimensionless  with 

-3  _4 

Pga   (i  =  1,2,3)  or  the  moment  made  dimensionless  with  p ga 

(i  =  4,5,6).   The  frequency  of  the  wave  excitation  is  denoted 

by  o,  and  6.  denotes  the  phase  shift  angle  of  the  i-th  component 

of  excitation  force  or  moment  in  relationship  to  the  incident 

wave.   (All  phase  shift  angles  are  measured  as  lag  positive 

and  in  relation  to  the  time  that  the  crest  of  the  undisturbed 

incident  wave  is  at  the  coordinate  origin.) 
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The  magnitude  of  the  complex  excitation  force  or  moment 
coefficient  occurring  in  (6.6)  is  defined,  as  in  Section  4, 


as 


1 (max) 

-3  0 

pga  n 


1=1,2,3;   C 


-   1(;"?  .1=4,5,6 

pga  n° 

* 


(6.7) 


where  n °  =  H/2a   denotes  the  dimens ionless  wave  amplitude,  H 

■k 

being  the  wave  height  and  a  the  characteristic  body  dimension. 

The  amplitude  of  the  excitation  forces  and  moments  are  denoted 

by  F.,    N  where  i  =  1,2,3  refers  to  the  three  components  of 
J        l (max)  r 

force,  and  i  =  4,5,6  refers  to  the  moment. 

As  the  body  responds  to  the  wave  excitation,  dynamic 
pressures  arise  due  to  the  motion  which  may  be  resolved  into 
two  components  of  force  (or  moment),  one  component  in  phase  and 
proportional  to  the  acceleration  of  the  body  and  a  second  in 
phase  and  proportional  to  the  velocity.   Equation  (4.10)  defines 
this  dimens ionless  force  (or  moment)  which  is  due  to  the  motion 
of  the  body.   The  dimens ionless  parameter  c. .(t)  is  defined 
in  (4.10)  as  the  force  (i=l,2,3)  or  moment  (i=4,5,6)  which 
is  caused  by  the  j-th  component  of  motion.   The  dimens ionless 
parameters  M. .  and  N. .  denote  the  added  mass  and  damping  coeffic 
ients  which  are  calculated  by  use  of  (4.6)  and  (4.7). 

The  final  contribution  to  the  force  (or  moment)  resulting 
from  the  surrounding  fluid  comes  from  the  hydrostatic  pressure. 
As  the  body  is  displaced  from  its  equilibrium  position,  forces 
and  moments  arise  which  are  proportional  to  the  body  displace- 
men  t  . 
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The  hydrostatic  pressure  increases  with  depth  according 
to  P  =  -P gy  and,  consequently,  the  i-th  component  of  hydro- 
static force  or  moment  resulting  from  this  pressure  variation 
and  acting  on  the  body  is  given  by  the  integral 


H 


F.   =  Pg 

1 


'a  x£y8idS' 


i  =  1,2,3; 


r,      H  _4 

F.   =  P  ga 

1 


y^Tyg.dS,  1=4,5,6    (6.8) 


where  y  =  y/a   denotes  the  dimens ionless  y  coordinate  of  a 

_   2 

point  on  the  immersed  surface,  and  dS =dS /a  ,  dS  being  an 

elemental  surface  area.   The  functions  g.  occurring  in  Eq .  (6.8) 
are  defined  as  follows: 


1     x 


:4  =  y  nz_Z  V 


h    '    ny 


g   =  z  n  -x  n  , 
s        x     z 


83  =  nz 


(6.9) 


g,  =  x  n  -y  n 
°6        y     x 


where  the  dimens ionless  body  coordinates  are  defined  as 

x'  =  x'/a,       y'  =  y'/a,       z'  =  z'/a,        d  =  d/a 

and  the  unit  normal  vector  on  the  immersed  surface  is  defined 

as   n(x',y',z';  Xn  ,  X0  ,  XQ  ,  X  .  ,  X  _  ,  X ,  )  =  in   +  jfn   +  kn  .   Eq.  (6.9)  is 

1   2.        3        4   D        b  xy       z 

the  same  definition  of  g.  given  in  (2.21f),  the  only  difference 

being  that  (2.21f)  is  written  in  terms  of  the  fixed  (x,y,z) 

coordinate  system  while  (6.9)  is  expressed  in  terms  of  the  body 

axes  . 

The  following  linearized  relationships  exist  between  the 

coorainates  of  the  fixed  reference  frame  and  the  body  coordinates 

for  a  given  point  on  the  immersed  surface  at  x',y',z': 

x  =  x'  +  Xcz  '  -  X,  y  '  +  X, 
->        6        1 


y  =  -d  +  y'  +  X,x'  -  X.z  *  +  X2 


z  =  z*  +  X.y '  -  Xcx'  +  X0 
4y      5       3 


(6.10) 
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Equation  (6.7)  represents  the  hydrostatic  force  (or  moment) 
acting  on  the  immersed  surface  including  the  buoyant  force 
which  must,  of  course,  just  balance  the  weight.   It  is  desired, 
however,  to  determine  the  i-th  component  of  force  associated 
with  a  j-th  (j  =  1,2, ...6)  component  of  displacement  of  the 
body,  j  =  1,2,3  denoting  linear  displacement  in  the  x',y',z' 
directions,  respectively,  and  j  =  4,5,6  denoting  angular  dis- 
placement about  x',y',z'  axes,  respectively.   For  the  linearized 
problem  this  force  (or  moment)  may  be  written,  using  Eq .  (6.8), 


as 


F.  . 


H 


3  F     ^ 

dr  xj  =  (1  or  ~a)  ?*l\  k.ff™±  ds  (6-11} 


j  --s 


(The  factor  1.0  in  brackets  in  (6.11)  is  applied  when  i=l,2,3  and 

a   is  appropriate  when  F . ,  denotes  a  moment  (i.e.),  i=4 , 5 , 6 . ) 
Defining,  further,  the  dimens ionless  parameter  K. .  as 


K 


« •  -// 


di(ygi)ds 


(6.12) 


the  dimens ionles s  force  (or  moment)  coefficient  k. .(t)  may  be 


written,  in  view  of  (6.11),  as 


k . . (t)  =  -K. .  X . (t)      (no  sum) 


(6.13) 


Equation  (6.13)  represents  a  form  appropriate  to  Eq.  (6.5)  for 
the  dimens ionles s  force  in  the  i-th  direction  caused  by  a  dis- 
placement in  the  j-th  degree  of  freedom. 
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The  following  expressions  are  obtained  for  the  hydrostatic 


force  coefficients  K. . 


K 


22 


=  -fb 


d  S  =  A  /a 
y         w 


K24  "  K42 


-/i-. 


dS 


K 


K 


,2. 
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(y'z'n   -  z  '  n  )dS 

J  z         y 


=  K64  =77sx'z'nydS 


K26  "  K62 


66 


•St 


,2 


(x'y'n   -  x'   n  )  dS 


(6  .14a) 
(6  .14b) 
(6  .14c) 
(6  .14d) 
(6  .14e) 
(6  .14f) 


in  which  A   denotes  the  waterplane  area  in  the  equilibrium 
w 

position.   For  floating  bodies  having  symmetry  with  respect  to 

the  x'-y'  and  y'-z'  planes  the  coefficients,  K0/  =  K. _  =  K.,  =  K, . 
J  j  t-  24     424664 

K26  =  K62  =  °* 

It  is  noted  here  that  the  expressions,  Eqs .  (6.14),  for  the 

hydrostatic  force  coefficients  are  in  rather  convenient  forms 
for  evaluation  as  a  part  of  the  same  numerical  scheme  discussed 
in  Chapter  5  which  involves  representation  of  the  immersed  sur- 
face by  a  large  number  of  nodal  points.   At  each  nodal  point  it 
is  necessary  to  specify  the  three  coordinates  of  its  location, 
the  three  components  of  the  unit  vector  normal  to  the  surface, 
and  the  elemental  area  of  the  facet.   This  information  is,  there- 
fore, available  in  convenient   forms  for  use  in  the  evaluation 
of  the  integrals  in  Eq.  (6.14)  in  the  computer  program. 
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If  the  mooring  lines  can  be  approximated  by  linear  elastic 

constraints  then  the  force  or  moment  coefficient  may  be  expressed 

in  terms  of  a  spring  constant  K' . .  and  displacement  X.  in  a  form 

13  J 


similar  to  Eq .  (6.12)  as 


k'  .  .(t)  =  -K»  .  .X.(t) 
ij  13  J 


( no  s  urn) 


(6.15) 


whe  re 


k-..(t)  - 


i  =  1,2,3;   j  =  1,2,  ...6:  force  in  i-di rec t ion/ p ga 


_3 


i  =  4,5,6;   j  =  1,2, ...6:  moment  in  i-di rec t ion/ pga 


_4 


V- 


and  K'  .  .  denotes  a  spring  constant  matrix  which  depends  on  the 

mooring  line  configuration,  stiffness  and  tension.   K' . .  is 

ij 

de  f  ined  as 


f9F 


K'  .  .  =  - 


K' . .  =  - 


3X 


9X. 


/  pga  , 


1=1,2,3 

j=l,2,. 


,   _4    1=4,5,6 
'PS*  ■   J-1.2....6 


(6.16a) 
(6.16b) 


in  which  F.  =  the  force  in  the  positive  i-th  direction  caused  by 

the  mooring  lines  when  i=i,2,3  and  F' .  =  the  moment  about  the 

i-th  body  axis  due  to  the  mooring  lines  when  1=4,5,6.   The 

parameters  X.  denote  dimens ionless  displacements  as  defined  by 
J 

Eq.  (6.3).   That  is,  when  i=l,2,3>X.  denotes  dimensionless  linear 


displacemen  ts 


X.  =  X. /a,   i  =  l,  2, 3 
11 


and  when  i=4,5,6  angular  displacements, 


X.  =  0 . ,     i=4, 5 ,6 
1     1 
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The  equations  of  motion  for  a  free  floating  body  may  now 
be  obtained  by  use  of  Eqs.  (6.4),  (6.5),  (6.6),  (4.10),  and 
(6.13)  as  follows : 


a  (m.  .  +  M..)X.  +  N..—  X.  +  (K..+K'..)Xj 
~    ij     iJ  3  13  g   J      iJ    iJ   J 


R  [  C,  e 


i6 


-lot 


(6.17) 
Furthermore,  the  body  response  in  the  j-th  degree  of  freedom  may 
be  expressed  in  the  form 


X.(t)  =  X!  R  [e^j    e-iat] 
3  J   e 


which,  when  substituted  into  Eq .  (6.17),  yields  the  complex 
equations  of  motion, 


(6.18) 


[-(m.  .+M.  .  )-iN  .  ,  +  (K.  .+K»  .  .) /v] 


X  .         I  C  .  I  .  x 
1  ±i>  a             116 
— 7T    e   J  =  ' '  e 


(6.19) 


where  the  frequency  parameter  v=  a  a/g. 

Equation  (6.19)  represents  six  equations  corresponding  to 
i=l,2,3...6.   The  repeated  j  index  denotes,  as  usual,  the  summa- 
tion over  the  six  degrees  of  freedom.   The  amplitude  ratio 

X°./r\°    denotes  the  ratio  of  the  amplitude  of  the  motion  to  the 

3       * 

amplitude  of  the  incident  wave  and  i/>.  denotes  the  phase  angle  of 
the  motion  (displacement)  in  relation  to  the  reference  condition 
of  the  crest  of  the  incident  wave  being  located  at  the  coordinate 
origin . 

It  will  be  recalled  that  for  bodies  possessing  x'-y'  and 
y'-z'  plane  symmetry,  m..=0  for  i#j  in  Eq.  (6.17).   Also,  M.. 
and  N. .  denote  the  thirty-six  element  added  mass  and  damping 
tensors . 
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It  can  be  shown,  however,  that 


K.  .  =K..,   M.  .  =M..,   K..  =  K.  . 
ij      Ji      iJ      ]i      iJ      ]i 


(6.20) 


so  that  only  21  different  values  exist  in  the  general  six  degrees 
of  freedom  problem. 

To  further  simplify  the  notation,  we  may  define  the  complex 
ma  tr ix , 

A..  =  -(m.  .+M.  .)  +  (K.  .+K'  .  .)/v-iN.  .  (6.21) 


the  complex  response  vector, 
.0 


Y 


X 


ii|< 


n' 


(6.22  ) 


and,     in    addition,     the    vector, 


C4  I  16  . 

i  1 

e 


v 


Then,  Eq.  (6.19)may  be  written  very  simply  as 

A    y,  =  B    i,j  =  1,2,  .  .  .6 
XJ       J     x » 


(6.23  ) 


Equation  ( 6 .  24  )  may  be  solved  by  matrix  inversion.   Once  the 

solution  to  Eq.  (6.21)  for  y.  is  obtained,  the  problem  is 

solved.   The  amplitude  and  phase  angle  of  the  response  can 

be  extracted  from  the  definition  of  y       given  in  Eq.  (6.22;. 
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7 .   HASKINDS  RELATIONS  AND  ENERGY  BALANCE 

Even  though  it  may  be  supposed  that  the  numerical 
solution  proposed  here  will  converge  upon  increasing  the 
number  of  partitions,  it  is  important  to  keep  the  partition 
size  large  (and  number  of  partitions  as  small)  as  accuracy 
considerations  will  permit  in  order  to  reduce  computer  time 
and  storage  requirements.   It  is,  for  this  reason,  impor- 
tant to  determine  the  effect  of  the  partition  size  on  accuracy 
so  that  practical  limits  may  be  established. 

One  method  of  verifying  the  accuracy  is  to  compare  the 
numerical  results  with  analytical  results  where  closed  form 
solutions  exist.   Although  valid,  this  approach  is  limited 
to  a  few  simple  shapes;  for  more  general  shapes  no  such  check, 
of  course,  exists. 

A  second  method  of  checking  the  validity  of  the  numerical 
results  involves  the  use  of  an  energy  balance  as  well  as  the 
use  of  the  so-called  Haskind's  relations.   Conservation  of 
energy  requires  that  a  balance  must  exist  between  the  energy 
required  to  oscillate  the  object,  and  the  wave  energy  trans- 
mitted across  some  control  volume  surrounding  the  object  but 
at  a  large  radial  distance.   Using  the  asymptotic  form  of 
the  Green's  function  given  in  (3.4a)  for  large  values  of  r, 
along  with  (3.1)  the  following  relationship  for  the  damping 
coefficient  is  so  obtained. 
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Nii  =  h         h(a'->)!v  £^ffn    V^I.OcoBhlaCh+TO] 

cos  (3-6) 


e"iapl 


,  2 
dS   d 


(7.1) 


where  p   =  /£ 2  +  n. 2 


-1 


and  6  =  tan  (  z,  /  £  )  .  This  relationship 
expresses  the  damping  coefficient  in  terms  of  the  far  field 
behavior  of  the  solution. 

A  relationship  somewhat  similar  to  (7.1)  known  as  Haskind's 
relations,  may  be  obtained  for  the  wave  force  (or  moment) 
coefficient.   That  is,  the  i-th  component  wave  force  (or 
moment)  coefficient  is  related  to  the  waves  produced  at 
infinity  by  the  body  oscillating  in  the  i-th  mode  such  that 


Ci    '   T^OThl     ffs    f^E.n.Occsh  [a(h+Il)]    e-laplc° 


s(3-tt-iJ^ 


dS 


(7.2) 


where  \\>     denotes  the  incidence  angle  as  indicated  in  Figure  1. 
A  form  of  this  relationship  between  the  wave  force  and  the 
waves  produced  at  infinity  by  the  same  body  oscillating  in 
otherwise  still  water  was  first  derived  by  Haskind  (14)  and 
later  reiterated  and  discussed  by  Newman  (3).   Equation  (7.2) 
may  be  considered  to  represent  a  form  of  the  Haskind's 
relations  as  extended  to  the  finite  depth  case.   The  details 
of  this  extension  which  involves  integration  by  method  of 
stationary  phase  are  given  by  Rao  and  Garrison  (15). 
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Equations  (7.1)  and  (7.2)  represent  relationships  for  the 

damping  and  wave  force  (or  moment)  coefficient  based  on  the 

behavior  of  the  far  field  solution.   A  comparison  of  these 

results  with  N..  and  C.  obtained  from  an  integration  of  the 
11       l 

pressure  over  the  immersed  surface,  i.e.,  as  obtained  from 
(4.2)  and  (4.7),  provides  a  convenient  and  valuable  self- 
check  on  the  accuracy  of  numerical  results.   These  results 
are  not  limited  to  special  configurations  and  may  be  applied 
to  arbitrary  shapes.   Equation  (7.2)  is,  however,  limited  to 
symmetry  with  respect  to  the  x-y  plane,  a  condition  which  is 
satisfied  by  most  practical  shapes. 

It  may  be  noted,  moreover,  that  for  the  special  case  of 
the  vertical  force  on  an  axisymmetric  body  a  very  simple 
relationship  between  C 9  and  N?9  may  be  obtained  from  Eqs.  (7.1) 
and  (7.2).   That  is,  in  view  of  the  axisymmetry  the  integrand 
in  Eq.  (7.1)  must  be  independent  of  6.   Accordingly,  the 
integration  with  6  may  be  accomplished  by  evaluating  the 
integrands  at  any  value  of  9,  say  6  =  tt   and  multiplying  by  2  tt  . 
The  right  hand  sides  of  Eqs.  (7.1)  and  (7.2)  then  become 
similar  and  the  following  relationship  may  be  derived  between 
the  wave  force  and  damping  coefficient  in  heave: 


N 


a.   sinh  (2  ah) 

22    2   2ah  +  sinh  (2ah) 


(7  .3) 
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For  the  case  of  infinite  depth,  Eq.  (7 .3)  reduces  to  the 
special  case  presented  by  Newman  (3) 


/2N~ 

v- 


C  i  =  /ZN22 
2    J— (7.4) 


The  numerical  evaluation  of  Eqs.  (7 .1)  and  (7  .2)  can  be 
easily  carried  out  using  the  same  numerical  scheme  outlined 
in  Chapter  5.   The  numerical  integrations  over  the  immersed 
surfaces  are  converted  to  summations  and  evaluated  using  the 
source  strength  function  obtained  from  solution  of  (5.2). 

8 .    COMPUTER  PROGRAM 

The  theoretical  development  outlined  in  Chapters  1-7  has 
been  coded  for  digital  computer  calculations  and  a  description 
of  that  program  is  presented  in  this  chapter.   The  general 
arrangement  of  the  numerical  procedure,  a  computational  flow- 
chart and  the  input-output  format  are  discussed. 

The  basic  requirements  established  for  the  program  were 
considered  to  be  as  follows: 

(a)  The  program  must  work  for  floating  bodies  of 
completely  arbitrary  shape,  but  when  either  one 
or  two  vertical  planes  of  symmetry  exist  the  pro- 
gram should  take  advantage  of  this  so  as  to  reduce 
the  input  data  and  calculations. 

(b)  Run  time  must  be  kept  as  small  as  possible  but 
without  losing  accuracy. 
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(c)   Computer  core  requirements  must  be  kept  to  a 
reasonable  value  even  for  very  complex  shapes 
where  a  large  number  of  panels  are  required. 
Flow  Chart  ' 

The  flow  chart  for  the  computer  program  is  given  in  Table  1 
It  should  be  noted  that  the  total  program  is  divided  into  two 
parts,  a  very  small  MAIN  program  which  is  used  only  for  dimen- 
sioning, inputing  parameter  and  calling  the  first  subroutine 
DYNRE3  •  This  arrangement  allows  all  of  the  subroutines  which 
have  variable  dimensioning  to  be  placed  on  files  if  desired 
and  only  the  small  MAIN  program  need  be  altered. 

The  subroutine OTECPR  controls  the  flow  of  the  calculations 
and  calls  essentially  all  of  the  other  subroutines.   It  first 
calls  the  subroutine  GDATA  which  in  turn  calls  the  subroutine 
CONFIG.   This  subroutine  reads  data  cards  describing  the 
immersed  surface  and  computes  the  unit  normals,  areas  and 
centroid  locations  of  each  panel.   GDATA  then  calls  the  sub- 
routine PLT  which  plots  the  panels  on  a  calcomp  plotter.   The 
CPU  time  requirement  to  this  point  in  the  program  is  small, 
and  therefore  the  program  could  be  terminated  at  this  point 
so  that  the  plot  output  could  be  checked.   Once  it  is  deter- 
mined that  the  input  data  is  correct  a  complete  run  could 
then  be  made. 

Next  the  subroutine  WAVLIN  is  called  by  OTECPR  .  The 
purpose  of  this  subroutine  is  to  compute  the  value  of  "a" 
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TABLE  I 
COMPUTER  PROGRAM  FLOW  CHART 


0TECPR1 


OTECPRIM 


Dimensions 
Input  parameters 


OTECPR 

Head 
mass  and 
s  n  r  3  n  e 
constant 
matrix 

caras 


PLT 


GDATA 


CONFIG 


TVAVLIN 


GCOEF 


OCINV 


SCAT 


OCINV2    -l 


Plots  the  immersed  surface  on  a 
Calcomp  plotter.  Input:  ALF  and 
BET,  the  azimuth  and  elevation 
angle  of  observer  and  R  the  dis- 
tance of  the  observer. 


Reads  data  cards  describing 
the  coordinates  of  the  corners 
of  the  nanels  on  the  immersed 
surface.  Computes  the  components 
of  the  unit  normal  vectors, 
areas  of  the  panels  and  coordin- 
ates of  the  centroids  of  the 
nanels . 


Computes  the  value  of  a  using 
the  "/ater  denth  and  wave  neriod 
in  Eq. (3.3f ) . 


Computes  parameters  used  in 

GREENS  including  the  trial 

and  error  solution  of  Eq.(3.4b) 


Inverts  ^--matrix  out  of  core. 


■^i*  (when  NSS=2) 

Inverts  (\-  matrix  out  of  core 
(when  NSS=1) 

Computes  Q( ,  /S  ,(Xx,C{u,(Yz        matrices 
and  stores  them  on  direct  access 
files  or  in  core  depending  on  the 
value  of  the  parameter,  NS . 


GREENS    1 


GREEN 


or 


GS I NG 


AVEVAL 
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Green's  function:  series  form. 
Green's  function:  integral  form. 

Comnutes  the  1/R  term  in  G. 
See  Appendix  B. 

ronnnips  c.  :is  riven  bv  Fm.C2.21  f) 


GCOEF 


OCINV2 


VELNOR 


SOURCE 


COMAT 


BASPR 


BAcCOE 


HASK 


REGPON 
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GREEN?  including  the  trial 
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ComDUtes  the  1/R  term  in  G. 

See  Appendix  B. 


Commutes  g.  as  given  by  Eq.(2.21f) 


Out  of  core  solution  of  Eq.(5.2) 

In  core  solution  of  Eq . ( 5 . 2 ) 

Computes  the  Dotential  at  the 
node  points  and  prints  the 
pressure  and  velocity  dist. 


Using;  the  oressure  distribution 
BASCOE  computes  the  excitation 
forces  and  moments,  added  mass 
and  damning  coefficients. 

Commutes  the  value  of  the  excita- 
tion forces  and  moments  and  added 
mass  and  damping  coefficients  by 
use  of  Haskind's  relations  and 
energy  balance. 

Solves  the  equations  of  motion  of 
the  floating  body, (i.e.),  solves 
Eq.  (  6.24)  . 
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Loops  through  :,N*VAVES"  times.  Each  time  through  the  loon  a  new 
wave  incidence  angle  is  read  at  the  same  oeriod  and  denth. 
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by  trial  and  error  solution  of  the  equation: 

v  =  a  tanh  (ah) 

The  subroutine  GCOEF  is  next  called  by  OTECPR.    The  func- 
tion of  this  subroutine  is  to  compute  certain  functions  which 
are  used  repeatedly  in  the  series  form  of  the  Green's  function. 
The  first  fifty  values  of  p   given  by  (3.4b)  are  computed  and 
these  values  are  used  to  compute  the  first  fifty  values  of 
cos  [  m,  (h+y ) ]  and  sin  [y,  (h+y)]  for  the  various  different  values 
of  y  corresponding  to  the  node  points  at  the  centroids  of  the 
panels  on  the  caisson.   The  values  of  these  functions  are 
stored  and  are  used  by  the  subroutine  GREENS  for  use  in  eval- 
uating the  series  form  of  the  Green's  function. 

The  subroutine  SCAT  is  next  called  by   OTECPR.   This  sub- 
routine evaluates  the  elements  of  the  a  and  B  matrices  and 
stores  them  on  FILE  #10  and  14,  respectively.   SCAT  calls  the 
subroutine  GREEN  and  GSING  or  GREENS  to  evaluate  these  matrices 
GSING  calls  AVEVAL  to  evaluate  the  1/R  and  3(l/R)/3n  term  in 
the  integral  form  of  the  Green's  function. 

The  program  has,  in  general,  two  options  for  solving  the 
matrix  equation  (5.2)  either  by  inversion  of  (a-I)  through 
OCINV  or  0CINV2  out  of  core,  or  by  solution  through  COMAT  in 
core.   When  the  inversion  out  of  core  option  is  selected, 
SCAT  calls  OCINV  or  0CINV2  to  invert  (a-I)  stores  (a-I)~   in 
the  file  space  originally  occupied  by  (a-I). 
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In  addition  to  the  calculation  of  a  and  g  matrices  and  the 

storing  of  these  matrices  out  of  core,  SCAT  also  generates  the 

derivatives  of  a  (i.e.,  the  matrices  a  ,  a   and  a   and  stores 

x    y       z 

them  on  Files  11,  12,  and  13). 

At  this  juncture  in  the  flow  of  calculations OTECPR   reads 
the  first  of  a  series  of  data  cards  which  specifies  the  incidence 
angle,  iJj  ,  of  the  incident  wave  relative  to  the  body  coordinates. 
All  of  the  computations  made  heretofore  are  dependent  on  the 
period  and  wave  height  but  are  independent  of  the  incidence  angle 
Thus,  computations  may  be  made  for  a  series  of  incidence  angles 
at  little  additional  expenditure  of  CPU  time.   The  remainder  of 
the  program  is,  therefore,  "looped"  NWAVES  times  with  different 
values  for  the  incidence  angle. 

The  next  subroutine  called  by  OTECPR  is  VELNOR.   This  sub- 
routine computes  the  velocity  induced  normal  to  the  caisson  by 
the  incident  wave.   More  specifically,  VELNOR  computes  the  vector 
given  by  (2.21f).   This  vector  represents  the  right  hand  side 
of  Equation  (5.2)  so  that  the  equation  can  be  solved  for  the 
source  strength. 

Depending  on  the  particular  option  selected  for  solving 
(5.2)  either  SOURCE  or  COMAT  is  called  next  in  OTECPR.   The 
function  of  SOURCE  is  to  multiply  (a-I)    on   g   in  order  to 
evaluate  the  source  strength  through  (5.2).   In  the  case  that 
(a-I)    has  not  been  determined  through  OCINV  or  0CINV2,  the 
subroutine  COMAT  is  called  which  solves  (5.2)  by  elimination. 
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The  subroutine  BASPR  is  called  next  by  OTECPR.   The  purpose 

of  this  subroutine  is  to  compute  the  value  of  the  potentials, 

(u_  +  u_) ,  u.,  u  ,  u  ,  u  ,  u   and  u,  at  the  node  points  on  the 
0712345        6 

immersed  surface.   Using  these  values  of  the  potentials,  BASPR 

computes  the  pressure  distribution  on  the  immersed  surface. 

This  subroutine  also  reads  the  FILES  #11,  12,  and  13  which  store 

a  ,  a   and  a   and  with  this  information  computes  the  three 
x    y        z 

components  of  velocities  in  body  coordinates  at  the  node  points 
on  the  immersed  surface.   It  should  be  noted  that  the  velocities 
are  used  only  for  a  check  on  the  solution  and  are  not  actually 
needed  in  the  computation  of  the  loads  and  response  of  the  floating 
body.   The  computation  of  velocity  can  be  surpressed  by  simply 
setting  the  parameter  NBA  =  1  in  the  MAIN  program;  when  NBA  =  2 
the  velocities  are  computed.   When  NBA  =  1  the  FILES  #11,  12  and 
13  are  not  used  at  all. 

The  subroutine  BASCOE  is  next  called  by  OTECPR.   The  function 
of  this  subroutine  is  to  re-compute  the  excitation  forces  and 
moments  by  use  of  the  Haskind's  relations  (Eq.   7.2)  and  re-compute 
the  damping  coefficients,  N..,  by  use  of  an  energy  balance  (Eq.  7.1) 
These  results  which  are  computed  on  the  basis  of  the  far-field 
potential  are  compared  with  the  same  results  computed  by  use  of 
the  surface  (near-field)  pressure  distribution  as  a  check  on  the 
validity  of  the  solution. 

Finally  OTECPR  calls  the  subroutine  RESPON.   This  subroutine 
takes  as  inputs  the  excitation  forces  and  moments,  added  mass  and 
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damping  coefficients,  mass  and  spring  constant  information  and 
then  solves  Equation  (6.24)  for  the  dynamic  response  of  the  body. 
Input  Data 

The  input  data  for  the  program  describes  the  geometry  of 
the  immersed  surface  up  to  the  mean  waterline,  specifies  the 
water  depth  and  wave  conditions  as  well  as  the  mass  and  mooring 
line  spring  constant  parameters.   Some  of  the  data  is  input 
through  data  cards  and  some  through  constants  in  the  DATA 
statements  in  the  MAIN  program. 

The  36  values  of  mass  m  and  mass  moments  of  inertia  of  the 

floating  body,  I-,-,,  etc.,  are  input  through  six  data  cards, 

x  y 

each  card  containing  six  values. 

The  spring  constants  which  account  for  the  effect  of  the 
mooring  lines  are  next  input.   Here  again,  there  are,  in  general, 
36  values  which  are  input  through  six  data  cards,  each  card 
containing  six  values. 

The  immersed  surface  of  the  body  is  described  through  a 
set  of  data  cards  which  give  the  three  coordinates  of  the  corners 
of  the  panels.   These  coordinates  are  measured  with  respect  to 
a  special  coordinate  system  used  for  inputing  the  data  only. 
These  special  "input  axes"  allow  the  user  of  the  program  to 
measure  the  coordinates  of  the  panel  corners  in  any  convenient 
reference  frame  and  then  the  subroutines  GDATA  and  CONFIG  shift 
the  origin  of  this  special  input  coordinate  system  to  any  desired 
location  by  use  of  the  inputs  XREF,  YREF,  ZREF  and  rotates  it 
by  the  angle  ANGREF. 
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Another  set  of  data  cards  following  the  ones  specifying 
the  location  of  the  corner  points  is  used  to  specify  which  four 
corners  constitute  a  given  panel.   Each  data  card  in  this  set 
contains  five  integers,  the  first  represents  the  panel  number 
or  index  and  the  following  four  integers  specify  the  indices 
associated  with  the  corners. 

The  last  set  of  data  cards  specify  the  different  values 
of  the  incidence  angle  of  the  wave.   There  are  "NWAVES"  of 
these  cards . 
Coordinate  Systems 

It  is  important  to  understand  the  coordinate  systems  used 
in  order  to  input  the  program  properly  and,  therefore,  these 
are  defined  in  the  following.   It  is  convenient  to  define 
three  coordinate  systems: 

(x",  y",  z"):  Input  axes 

(x',  y',  z'):  Body  axes 

(x,  y,  z)     :  Inertial  axes 

The  input  axe s  are  used  to  input  the  location  of  the  corners 
of  the  panels.   This  set  of  axes  may  be  located  at  any  convenient 
point  and  orientation  with  respect  to  the  hull  but  the  y"-axis 
is  set  vertical  with  positive  y"  measured  upward.   The  second 
two  coordinate  systems  also  have  their  y-axes  pointing  upwards 
and  the  location  of  their  origins  and  orientation  is  specified 
by  the  five  parameters,  XREF,  YREF,  ZREF,  ANGREF  and  ECG. 
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The  inertial  axes  is  located  such  that  the  x'-z'  plane 
represents  the  mean  water  level.   It  is  located  at  the  point 
x"  =  XREF,  y"  =  YREF  and  z"  =  ZREF  in  the  input  reference  frame 
such  that  the  center  of  gravity  of  the  body  is  on  a  vertical 
line  passing  through  the  origin.   The  body  axes  is  considered  to 
be  aligned  with  the  inertial  axes  when  the  body  is  in  its  equili- 
brium position  except  that  it  is  shifted  to  the  point  (0,  ECG,  0) 
in  the  (x,y,z)  reference  frame.   When  the  center  of  gravity  lies 
above  the  mean  waterline  ECG  represents  a  positive  number  and 
when  the  center  of  gravity  lies  below  the  mean  water  level  ECG 
will  be  negative. 

Figure  4  shows  the  three  coordinate  systems  and  their 
relationship  one  to  another.   The  reason  for  the  introduction 
of  the  input  coordinate  system  is  strictly  for  convenience  in 
inputing  data.   In  the  case  of  certain  configuration  it  may  be 
convenient  to  measure  the  location  of  the  corners  of  the  panels 
in  a  coordinate  system  attached  to  the  hull  at  a  reference 
point  which  may  differ  from  either  the  mean  waterline  or  the 
center  of  gravity. 

Positive  values  of  XREF,  YREF,  ZREF  and  ANGREF  are  shown  in 
Figure  4.   The  wave  incidence  angle   ^  }  is  also  shown  in  the 
figure  and  is  defined  relative  to  the  body  and  inertial  axes 
as  indicated . 
Corner  Node  Point  Inputs  -  No  Symmetry: 

It  would  be  unusual  that  a  floating  structure  would  possess 
no  symmetry  whatsoever  but  the  program  has  the  capability  of 
dealing  with  such  general  cases.   When  no  symmetry  is  present 
the  complete  immersed  surface  of  the  hull  must  be  covered  with 
panels . 
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FIGURE  4   DEFINITION  OF  COORDINATE  SYSTEMS 


Corner  Node  Point  Input  -  One  Plane  of  Symmetry: 

When  the  immersed  surface  has  one  vertical  plane  of 
symmetry  and  the  center  of  gravity  lies  in  this  plane,  then  only 
half  of  surface  must  be  represented  by  panels.   In  such  a  case 
the  input  axes  may  be  placed  in  any  desired  manner  relative  to 
the  immersed  surface  but  the  parameters  XREF,  YREF,  ZREF  and 
ANGREF  must  be  selected  such  that  the  plane  of  symmetry  is  coin- 
cident with  the  x1  -  y'  and  x  -  y  planes  and  such  that  0(x',y',z') 
is  coincident  with  the  center  of  gravity.   The  single  symmetry 
case  is  depicted  in  Figure  5. 
Corner  Node  Point  Inputs  -  Double  Symmetry: 

If  the  hull  possesses  two  planes  of  symmetry  it  is  necessary 
to  input  data  cards  specifying  the  panel  corners  on  one  quarter 
of  the  hull  only.   Here  again  the  input  coordinate  system  i s 
used  to  input  the  coordinates  of  the  corners  of  the  panels  and 
it  may  be  positioned  arbitrarily  with  respect  to  the  hull. 
However,  XREF,  YREF,  ZREF,  ANGREF  have  to  be  input  such  that 
0(x,y,z)  and  0(x',y',z')  lie  on  the  intersection  of  the  two 
planes  of  symmetry  and  ECG  must  be  such  that  0(x',y',z')  is 
coincident  with  the  center  of  gravity.   The  result  is  that  the 
(x,y,z)  and  (x',y',z')  are  placed  relative  to  the  floating  body 
as  shown  in  Figure  6.   The  portion  of  the  immersed  surface  input 
falls  in  the  positive  x-z  quadrant  in  the  body  and  inertial 
axes  as  indicated  in  the  figure. 
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FIGURE  5   SINGLE  SYMMETRY 


6  7 


y" 


Input  -^  Surface. 


xfief* 


V 


I    J/ 


& 


/ 


'II 


MGURE    6       DOUBLE    SYMMETRY 
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Although  in  the  above  two  cases  where  the  body  possesses 
one  or  two  planes  of  symmetry  the  input  axis  was  defined  such 
that  XREF,  YREF,  ZREF,  and  ANGREF  were  non-zero,  the  usual  case 
would  be  to  input  the  data  with  input  coordinate  axes  aligned 
with  the  body  and   inertial  axes.   This  is  generally  convenient. 
However,  it  is  not  usually  convenient  to  place  the  origin  of 
the  input  axes  at  either  the  mean  water  level  or  center  of 
gravity.   It  is  generally  more  convenient  to  set  the  origin  at 
the  level  of  the  bottom  if  the  hull  has  a  flat  bottom  so  that 
YREF  is  normally  non-zero. 
Input  Format 

The  inputs  to  the  computer  program  are  of  three  types: 
(1)  specifying  certain  constants,  (2)  specifying  the  dimensions 
of  the  arrays,  and  (3)  specifying  the  geometry  of  the  immersed 
surface.   These  inputs  and  their  format  are  discussed  in  the 
following : 
Input  Constants: 

The  constants  which  are  input  by  way  of  DATA  statement 
cards  are  defined  in  terms  of  both  the  metric  and  English  system 
in  the  following.   It  is  noted,  however,  that  either  one  system 
or  the  other  may  be  selected  but  no  mixing  of  systems  is  allowable. 

PER Wave  period  and  period  of  the  motion  of  the  body  (sec.  ) 

H Mean  water  depth  (ft.  or  meters) 

ABAR Reference  dimension  of  the  hull.   Any  value  will 

do  but  normally  the  radius  of  the  hull  or  the 
half-length  is  used  (ft.  or  meters). 
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3         3 
RHO Density  of  the  water  (slugs/ft   or  kg/in  ). 

2 
G Acceleration  of  gravity  (32.174  ft/sec   or 

9. 8m/sec2)  . 

XREF x"  coordinate  of  the  origin  of  the  body  axes 

(measured  in  the  input  axes )  (f  t .  or  meters). 
(See  Figure  4 ) . 

YREF y"  coordinate  of  the  origin  of  the  body  axes 

(measured  in  the  input  axes ) ( f t  or  meters) 
(See  Figure  4) . 

ZREF z"  coordinate  of  the  origin  of  the  body  axes 

(measured  in  the  input  axes ) ( f t  or  meters) 
(See  Figure  4) . 

ANGREF. .. Clockwise  rotation  angle  of  the  body  axes 

(as  viewed  from  above)  (Degrees )  See  Figure  4. 


ECG 


Elevation  of  the  center  of  gravity  of  the  body 
(y)  measured  in  the  inertial  axes  located  with 
x-z  plane  on  the  mean  waterline. 
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RMIN The  value  of  RMIN  is  used  for  two  different  tests 

of  the  distance  between  panels  i  and  j  in  computing 
the  values  of  a.,  and  8...   The  use  of  RMIN  as  well 
as  ARMIN  is  bestJ  demonstrated  by  the  following 
flow  chart: 
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The  value  of  RMIN  is  used  to  test  R  for  purposes 
of  selecting  GREEN  or  GREENS.   In  subroutine  GSING 
it  is  used  to  determine  whether  or  not  AVEVAL  will 
be  called  to  integrate  the  1/R  and  derivative  of 
1/R  in  the  integral  form  of  the  Green's  function. 

The  value  of  RMIN  should  be  kept  as  small  as 
possible  consistent  with  adequate  accuracy.   Figure  Bl 
in  Appendix  B  indicates  that  RMIN  >_2  /a~S   (where  AS 
represents  the  area  of  a  typical  panel)  represents 
a  reasonable  value  fo  RMIN.   (ft  or  meters). 


NC 


Integer  denoting  the  number  of  panel  corner  points 
on  the  portion  of  the  hull  which  is  input.   This 
integer  represents  the  actual  number  of  data  cards 
read  by  CONFIG  which  describe  the  location  of  the 
corner  nodes  on  the  hull. 


NB 


NBA, 


Integer  denoting  the  total  number  of  panels  on 
the  complete  hull.   The  panels  are  assigned  indices 
running  1  through  NB .   Figure  7  shows  the  numbering 
system  when  the  hull  has  either  single  or  double 
symmetry . 

Integer  having  the  value  of  either  1  or  2.   If 
NBA  =  1  the  velocity  components  at  the  node  points 
on  the  hull  are  no  t  computed  and  direct  access 
files  #11,  12  and  13  are  not  used.   If  NBA  =  2  the 
velocity  components  on  the  immersed  surface  are 
computed.   It  should  be  noted  that  the  velocities 
are  not  actually  used  in  the  calculations  of  the 
pressures,  excitation  forces  and  moments,  added 
mass  and  damping  coefficients  or  the  resulting 
dynamic  response. 
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FIGURE  7  NUMBERING  SYSTEM  FOR  THE  CASE  OF  SINGLE  AND  DOUBLE  SYMMETRY 
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NA Number  of  panels  defined  by  the  NC  corner  nodes 

on:    (a)  the  complete  hull  for  the  no  symmetry 
case,  (b)  half  of  the  hull  when  the  hull  has  one 
plane  of  symmetry,  (c)  one  quarter  of  the  hull 
when  the  hull  has  two  planes  of  symmetry.   In  the 
case  of  (a),  (b)  and  (c)  NA  will  have  the  values 
NB,  NB/2  and  NB/4,  respectively. 

NEQ The  (a-I)  matrix  as  indicated  in  Eq .  (5.2)  is  a 

square  matrix  of  size  NB  x  NB  in  all  cases  except 
when  NSS  =  2,  in  which  case  it  is  NB/2  x  NB/2. 

In  all  cases,  however,  the  matrix  is  placed 
on  direct  access  FILE  #14  and  is  brought  into  core 
one  piece  at  a  time.   The  size  of  the  block  that 
is  brought  in  is  always  NEQ  (number  of  equations) 
by  NB.   Thus,  it  is  most  efficient  to  make  NEQ  as 
large  as  possible.   However,  the  block  size  is 
limited  (about  7200  bytes  on  an  IBM  360)  so  NEQ 
may  have  to  be  small  when  NB  is  large  so  that  the 
block  size  remains  within  the  limits  for  the  direct 
access  device  used. 

In  general,  NEQ  may  have  any  value  equal  to 
or  greater  than  unity.   The  general  rule  is  to  make 
NEQ  as  large  as  possible  consistent  with  the  maximum 
block  size  allowable. 

NWAVES ...  Integer  denoting  the  number  of  wave  directions 

which  will  be  considered  in  the  run.   There  must 
be  NWAVES  data  cards  giving  the  incidence  angle 
(in  degrees)  of  the  particular  cases  to  be  calcu- 
lated. 

NSS Integer  having  either  the  value  of  1  or  2.   The 

normal  value  for  NSS  is  1  but  under  the  special 
case  where  the  wave  direction  is  aligned  with  the 
x'-y'  plane  (0  or  180  degrees)  and  this  plane 
represents  a  plane  of  symmetry,  then  NSS  may  be 
set  to  2.   The  program  will  then  take  advantage 
of  the  fact  that  the  source  strength  function  on 
the  positive  z'-half  of  the  hull  is  the  same  as 
at  the  corresponding  points  on  the  negative  z'- 
half.   Accordingly,  the  (a-I)  matrix  will  be  only 
(NB/2)  x  (NB/2)  rather  than  NB  x  NB  and,  therefore, 
a  savings  in  CPU  time  will  be  realized  during  inversion 

NS Integer  which  is  normally  set  to  1.   However, 

under  the  special  case  where  NSS  is  set  to  2  the 
(a-I)  matrix  may  be  small  enough  to  store  in  core. 
If  this  is  the  case,  NS  may  be  set  equal  to  the 
integer  NB/2.   The  (a-I)  matrix  will  then  be  stored  in 
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the  ALS  matrix  and  inverted  through  subroutine 
COMAT.   It  is  always  less  time  consuming  to  take 
the  NS  =  NB/2  option  when  possible  provided  the 
storage  requirements  do  not  exceed  that  available 
since  the  in-core  inversion  is  faster  than  the 
out-of-core  inversion. 

Direct  Access  Files: 

The  DEFINE  FILE  statements  in  the  MAIN  program  specify  the 

block  size  and  number  of  records  on  the  direct  access  files  //10-14 

as  f ol lows : 

DEFINE  FILE  10  (a,b,U,LK) 

DEFINE  FILE  11  (a,b,U,LK) 

DEFINE  FILE  12  (a,b,U,LK) 

DEFINE  FILE  13  (a,b,U,LK) 

DEFINE  FILE  14  (c,d,U,LK) 

where,  for  example,  in  #10,  a  =  number  of  records  and  b  =  number 

of  words  in  the  record.   The  values  of  a  and  b  must  be  input  for 

a  given  set  of  data  and  conditions  as  follows: 

Number  of  records  - 

a  =  NA 

c  =  NB/NEQ  (rounded  upwards  to  the  nearest  integer) 

Block  size  -  (in  number  of  words) 

b  =  2xNB 

d  =  2xNEQxNB 

For  the  special  case  when  NSS  =  2  : 

c  =  (NB/2)/NEQ  (rounded  upwards  to  the  nearest  integer) 

In  addition  to  the  DEFINE  FILE  statements  it  is  necessary 

to  specify  adequate  and  consistent  disk  space  on  JCL  input. 


The  factor  of  2  accounts  for  the 
fact  that  the  words  are  complex. 
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For  efficient  operation  it  is  important  that  the  five  files 

#10-14  be  placed  on  separate  units  when  all  are  being  used  because 

these  five  files  are  written  upon  simultaneously  in  the  same 

DO  LOOP. 

When  NBA  =  1  files  // 1 1 ,  12  and  13  are  not  used.   In  this 

special  case  it  is  possible  to  set  the  file  space  required  for 

#11-13  to  some  nominal  value. 

Summary  of  Program  Options 

In  view  of  the  fact  that  the  program  has  several  options 

depending  on  the  incidence  angle  and  symmetry  of  the  immersed 

surface,  these  are  summarized  in  the  following.   The  various 

options  are  completely  independent  of  the  values  of  NBA  and  NEQ. 

I.   Hull  has  no  planes  of  symmetry: 

NA  =  NB 

NS  =  1 

NSS  =  1 

The  (a-I)  matrix  stored  on  file  #14  will  be  NBxNB.   The 

matrix  will  be  inverted  by  use  of  an  elimination  scheme 

through  subroutine  OCINV.   This  option  represents  the 

most  CPU  time  consuming  method  and  would  be  used  only 

if  the  hull  had  no  planes  of  symmetry  with  respect  to 

the  body  axes.   The  wave  incidence  angle  may  be  arbitrary. 

II.   Hull  has  either  one  or  two  planes  of  symmetry: 

NA  =  NB/2  (input  data  for  +  z  '  half  of  the  immersed  surface 
when  hull  has  one  plane  of  symmetry.) 

NA  =  NB/4  (input  data  for  +  x  '  and  +z '  quarter  of  the  hull 
when  the  hull  has  two  planes  of  symmetry). 
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(a)  NS  =  1 
NSS  =  1 

The  program  will  take  advantage  of  the  symmetry  in 
computing  the  matrices  but  the  (a-I)  matrix  will 
be  NB  x  NB  and  will  be  inverted  in  OCINV.   The  wave 
incidence  angle  may  be  arbitrary. 

(b)  NS  =  1 
NSS  ■  2 

The  program  will  take  advantage  of  the  symmetry  in 
computing  the  large  matrices.   The  incidence  angle 
must  be  set  to  either  0  or  180  degrees  so  that  the 
wave  is  aligned  with  the  x'-y'  plane.   The  assumption 
is  made  that  for  this  case  the  source  strength  on 
the  -z'  -half  of  the  hull  is  the  same  as  on  the  +  z  ' 
-half.   Accordingly,  the  matrix  stored  on  FILE  #14 
will  be  NB/2  x  NB/2  and  will  be  inverted  out  of  core 
by  0CINV2. 

(c)  NS  =  NB/2 
NSS  =  2 

The  (a-I)  matrix  stored  on  FILE  #14  will  be  of  size 
NB/2  x  NB/2.   This  matrix  is  eventually  placed  in 
core  in  the  array  ALS(NS,NS)  which  for  this  option 
has  been  dimensioned  NS  x  NS,  i.e.,  the  size  of  the 
(a-I)  matrix.   The  fact  that  NS  is  other  than  unity 
signals  the  computer  to  solve  Eq.  (5.2)  in  core  by 
use  of  the  subroutine  COMAT.   For  this  option 
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(NSS  =  2)  the  incidence  angle,  \\>     ,  must  be  0  and/or 
180  degrees  so  that  the  wave  direction  is  aligned 
with  the  plane  of  symmetry. 
Dimensioning 

The  general  principle  to  be  followed  in  the  dimension 
statements  in  the  MAIN  program  is  to  make  the  dimensions  exact ly 
equal  to  the  size  of  the  array.   In  the  case  of  two-dimensional 
arrays  it  is  absolutely  necessary  that  they  be  dimensioned 
exactly  correctly.   The  dimensions  cannot  be  either  larger  or 
smaller  than  the  actual  array.   Even  in  the  case  of  several  of 
the  one-dimensional  arrays  this  is  essential  and,  therefore, 
as  a  general  principle  the  arrays  should  always  be  dimensioned 
to  exactly  the  size  of  the  array  as  given  in  the  following  with 
the  exception  of  SH,  CH,  SINAMU  and  COSAMU  arrays. 

The  MAIN  program  contains  all  of  the  dimension  statements 
which  must  be  altered.   It  is  unnecessary  to  make  any  changes 
in  the  subroutines  because  all  arrays  are  variable  dimensioned. 
The  arrays  in  the  MAIN  program  have  the  following  dimensions: 

ALPHA  (NEQ,NB) 

ALS  (NS,  NS) 

AV(NB),  BV(NB),  CV(NB),  DV(NB),  EV(NB) 

F(NB,7),  HH(NB,7),  U(NB,7) 

XB(NA),  YB(NA),  ZB(NA),  XNB(NA),  YNB(NA),  ZNB(NA),  DSB(NA) 

XC (NC) ,  YC(NC) ,  ZC(NC) 

XP(NC),  YP(NC),  ZP(NC) 

N0DM(NA, 4) 

IC0M(NS, 3) 

SH(a),  CH(a),  SINAMU(b),  COSAMU(b),  where 
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a  _>  number  of  different  levels  of  node  points  on  the 
immersed  surface.   That  is,  "a"  should  be  greater  than 
the  largest  number  printed  under  the  column  labeled  LEVEL 
NO.  on  the  output.   A  value  of,  say  20,  would  cover  essen- 
tially all  conceivable  cases. 

b  >    50  x  a  (for  example,  b  =  1000  if  a  =  20) 
(Note:   In  the  case  of  SH,  CH,  SINAMU  and  COSAMU  the 
dimensions  may  be  greater  than  or  equal  to  the  actual  array 
size) . 

Input  Data  Cards 

Mas  s  Mat  r ix : 

The  first  set  of  data  cards  are  used  to  input  the  mass 

matrix.   There  are  six  cards  in  this  set;  each  card  contains 

six  numbers.   The  format  for  each  card  is  6F10.1  and  the  units 

2 
are  slugs  and  slug-ft   when  the  English  system  is  used  or  kg 

2 
and  kg-m   when  metric  units  are  used. 

The  first  card  contains  the  values  of  m..  ..  ,  m,  9 ,  m..  _  ,  m  ,  , 

mn  c  and  mn  ,     the  second  card  contains  values  of  m„n  m_0,  ni„_,  m„. 
±_>        lb  ZL       Z  Z  Z  3  Zh 

m„,-,  m„,,  etc.  where 
25    26 

m1 ,  =  m    =  m„_  =  m  =  mass  of  the  structure  (slugs  or  kg) 

m12  =  m21  =  m13  =  m^    =    m23  =  m^    =    0 

2 
m,,  =  1-,-,=  moment  of„inertia  about  the  x'  axis  (slugs-ft 

o  r  kg  -  m  ) 

2 
.  =  moment  of  inertia  about  the  v'  axis  (.slugs-tt 

y  y 


mcc  =  I-,-,  =  moment  of  inertia  about  the  y'  axis  (slugs-ft 
55     y  y  2X 


or  kg-m  ) 
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m 
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=  I-,-,   =  moment  of  intertia  about  the  z'  axis  (slugs  -  ft 
or  kg  -  m^  ) 


(Note,  the  x',  y'  and  z'  axes  are  centroidal  axes  by  definition) 


m 


=  m, 


4  5     5  4 


=  -I-,-,  =  negative  product  of  inertia  (slugs  - 


x  y 


ft   or  kg  =  mz) 


m.   =  m,.  =  -I-,-,  =  negative  product  of  inertia  (slugs  - 

q-D        (11         XZ        r  *.  9         i  ^  \ 

ftz   or  kg  -  m  ) 


m  ,  =  m   -- I- , -  , 
56     65    y  z 


=  negative  product  of  inertia  (slugs  - 
ft   orkg-m) 


Spring  Constant  Matrix : 

The  second  set  of  data  cards  are  used  to  describe  the 
elastic  constraints  of  the  body  due  to  mooring  lines.   In  the 
most  general  case,  for  displacement  in  each  of  the  six  degrees 
of  freedom  a  force  (i  =  1,2,3)  and  moments  (i  =  4,5,6)  can  occur. 
Thus,  in  general,  there  can  be  36  elastic  coefficients  to  describe 
the  effect  of  the  mooring  lines. 

The  definition  of  the  elements  of  the  spring  cons  tan t 
ma  t  r i  x  is  as  follows: 


ij 


3f! 

1 
3X  . 


where 


,    Iforce  in  the  i-direction  (i=l,2,3)(lb  or  N) 
i     [Moment  about  i-axis  ( i=4 , 5 , 6 ) ( f t-lb  or  N-m) 

Linear  displacement  in  the  j -di r ec t ion ( j =1 , 2 , 3 ) 
(ft  or  meters ) 


Xj  -< 


Angular  displacement  about  j-axis  (j-4,5,6) 
(radians ) 


v 


Examp le : 

The  case  of  a  simple  spring  restraint  is  indicated  in  the 


following  sketch: 


WVWW p 

k, 


It  is  possible  to  derive  expressions  for  the  values  of  K ' . .  by 

starting  with  the  expressions  for  the  two  components  of  force 

and  moment: 

F^  =  -k1(X1-X1  ) 
o 


F»   =  -k2(x2-x2  +  a±x  ) 

o 


-k3(X2-X2   -  £2X6) 
o 


F'6  =  -k2£1(X2-X2   +  ^XX6) 

o 


+  k3  £2(X2-X2  -^2X6) 
o 


where  the  subscript  "o"  denotes  the  equilibrium  position  of  the 
body.   By  differentiation,  the  spring  constant  coefficients 

0 

K'  .  .  are  defined  as: 

in.   oo      K.  «  •  K.  — 

2  2 

K'66  =  k2£l   +  k3  £2 

£,61  =  ^16'  =  ° 
K'12  -  K'2l  =  0 

R,26  =  R,62  =  k2£l  "k3£2 
where  k1 ,  k„,  and  k„  denote  the  spring  constants  of  the  three 
springs  defined  in  the  usual  manner  in  terms  of  force/displacement 
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Hull  Geometry: 

The  next  set  of  data  cards  to  be  read  are  the  ones  which 
describe  the  immersed  surface  of  the  hull.   This  set  of  cards 
gives  the  coordinates  of  the  four  corners  of  the  panels  and  there 
will  be  NC  of  these  cards.   When  there  is  no  symmetry  in  the 
immersed  surface,  the  complete  hull  must  be  described;  when  the 
hull  has  a  single  plane  of  symmetry  (x'-y'  plane),  then  half  the 
surface  must  be  described,  and  when  the  immersed  surface  has  two 
planes  of  symmetry,  one  quarter  of  the  hull  must  be  described. 
When  inputing  half  the  immersed  surface,  the  x'-y'  plane  must 
represent  the  plane  of  symmetry  and  when  inputing  one  quarter  of 
the  hull,  the  x'-y'  plane  and  y'-z'  planes  must  represent  planes 
of  symmetry.   In  the  case  of  single  symmetry  the  +  z'  part  of  the 
immersed  surface  is  input, and  in  the  case  of  double  symmetry 
the  +x'  and  +z'  quadrant  is  input. 

An  example  of  a  subdivided  surface  with  two  planes  of  symmetry 
is  shown  in  Figure  8.   On  the  quarter  of  the  hull  to  be  described 
there  are  60  corner  node  points  denoted  and,  therefore,  for  this 
example,  NC  =  60.   It  should  be  noted  that  all  panels  must  have 
four  sides;  however,  the  length  of  the  sides  and  orientation  may 
be  arbitrary,  i.e.,  the  panels  need  not  be  rectangular  or  square. 

The  coordinates  of  the  NC  (60  in  the  example)  panel  corner 
points  must  be  input  by  use  of  NC  data  cards  punched  in  the  fol- 
lowing format : 

Cols:  1-10:   The  corner  node  index  (1-60  in  the  example)  (Integer, 
right  ad j  us  ted)  . 
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Cols:  11-20:  x"  coordinate  of  the  node  (ft  or  m) . 
Cols:  21-30:  y"  coordinate  of  the  node  (ft  or  m) . 
Cols:   31-40:  z"  coordinate  of  the  node  (ft  or  m) . 

Correspondence  Table: 

The  next  set  of  data  cards  will  define  the  sequence  of  four 
corner  nodes  which  make  up  a  given  panel.   There  are  five  integers 
on  each  of  these  data  cards,  the  first  represents  the  panel  index 
and  the  next  four  integers  correspond  to  the  four  corners  of  the 
given  pane  1 . 

There  will  be  NA  of  these  data  cards  where  NA  denotes  the 
number  of  panels  on  the  portion  of  the  immersed  surface  described 
by  data  cards.   NA  will  have  the  value: 

NA  =  NB  (no  symmetry) 

NA  =  NB/2  (single  symmetry) 

NA  =  NB/4  (double  symmetry) 
where  NB  denotes  the  total  number  of  panels  on  the  complete 
immersed  surface.   In  the  example  described  in  Figure  8: 

NC  =  60 

NB  =  188 

NA  =  47 
In  this  case  the  hull  has  double  symmetry  so  NA  =  NB/4  -  188/4  =  47 

Each  of  the  NA  cards  will  specify  which  of  the  four  points 
compose  a  given  panel  and  are  punched  according  to  the  following 
format  (all  right  adjusted  integers): 
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Cols:   1-10:  Panel  number  (1  through  47  in  the  example) 
Cols:   11-20:  First  corner  node  index 
Cols:   21-30:  Second  corner  node  index 
Cols:   31-40:  Third  corner  node  index 
Cols:   41-50:  Fourth  corner  node  index 

The  four  corner  nodes  are  specified  in  sequence  as  encountered 
in  "walking"  around  the  panel  on  the  water  side  in  such  a  way 
that  the  interior  of  the  panel  area  remains  on  the  right  hand  side. 
Another  equivalent  way  of  stating  this  is  to  number  in  sequence 
moving  around  the  panel  clockwise  as  viewed  from  the  water  side 
of  the  hull  surface.   The  starting  corner  point  is  arbitrary. 

Wave  Incidence  Angle  Data  Cards: 

The  final  set  of  data  cards  are  those  specifying  the  different 
wave  directions  to  be  considered  in  the  run.   Since  the  (a-I) 
matrix  needs  to  be  constructed  and  inverted  only  once,  it  requires 
little  additional  CPU  time  to  consider  several  values  of  incidence 
angle  in  addition  to  the  initial  one.   It  should  be  remembered, 
however,  that  when  NSS  =  2  the  wave  must  be  aligned  with  the  x' 
axis  so  the  incidence  angle  can  only  be  0  and/or  180  degrees  only. 

The  format  for  the  cards  are: 
Cols:   1-10:  Incidence  angle,  ip  ,  in  degrees. 

Computer  Print-Out 

The  computer  print-out  contains  some  of  the  input  data  as 
well  as  the  computed  results. 
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The  coordinates  of  the  corner  nodes  are  printed  as  read-in 
except  that  they  are  shifted  to  the  inertial  axes  (x,y,z)  with 
origin  at  the  mean  water  level  on  a  vertical  passing  through  the 
e.g.  of  the  body.   The  units  are  the  same  as  input. 

The  correspondence  table  is  also  printed  as  read  in  for 
purposes  of  double  checking. 

The  centroids  (node  points)  on  the  panels  are  computed  as 
well  as  the  components  of  the  unit  normal  vectors  and  areas  of 
the  NA  panels  input.   These  results  are  printed  under  the  heading, 
"COORDINATES  OF  CENTROIDS  OF  PANELS,  COMPONENTS  OF  N  AND  PANEL 
AREAS. " 

The  displacement  of  the  immersed  surface  is  computed  by 
three  different  methods: 

V  =  ff  xn  dS 

JJs      x 


¥  -ffsy   ° 


dS 

s*     y 


-lis- 


dS 
z 


The  results  for  the  volume  computed  in  the  three  different  manners 
represents  a  good  quick  check  on  the  accuracy  of  the  input  data. 
If  the  input  data  is  correct,  the  three  results  should  be  the  same 

In  a  similar  fashion  the  total  surface  area  as  well  as  the 
area  projected  in  both  the  x  and  z  directions  are  computed  as  a 
quick  check  on  the  input  data. 

Next  the  program  goes  through  the  NA  panel  node  points  and 
assigns  an  index  to  each  different  elevation.   All  of  the  node 
points  at  elevation  -10.0,  for  example,  may  be  assigned  the  index 
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1  and  all  at  -5.0  would  be  assigned  2,  etc.   The  panel  node  point 
index  and  its  corresponding  "level  number"  is,  therefore,  printed. 
Caution :   largest  number  printed  under  "LEVEL  NO."  should  not 
exceed  the  dimension  of  the  arrays  SH  and  CH .   The  maximum  value 
would  typically  be  3  to  10.   A  dimension  of  20  to  30  should 
include  all  cases. 

The  "WAVE  NUMBER"  is  next  printed.   This  is  simply  an  index 
which  indicates  that  a  new  wave  direction  is  being  considered. 
The  number  printed  here  will  range  from  1  through  "NWAVES." 

The  next  block  of  information  in  the  print-out  is  the 
f o llowing : 

a  =  2TTa/L  =  2 tt  .  hull  radius/wave  length. 

v  =  (2WT)2i/g 

h  =  water  depth  (ft  or  m) 

T  =  wave  period  (sec) 

ip  =  wave  incidence  angle  (degrees)  See  Figure  1 

NB  =  total  number  of  panels  on  total  hull  surface. 

The  next  block  of  data  printed  represents  the  pressure  dis- 
tribution associated  with  the  motion  of  the  body  in  its  six  degrees 
of  freedom:   surge  (1),  heave  (2),  sway  (3),  roll  (4),  yaw  (5), 
and  pitch  (6).   The  subscript  (7)  denotes  wave  interaction  with 
the  fixed  body . 

For  each  mode  (six  degrees  of  freedom  1-6  and   7   wave  inter- 
action with  the  fixed  hull)  the  dimensionless  pressure  amplitude 
coefficient  PI,  P2,  etc.  and  pressure  phase  angle  PHI,  PH2,  etc. 
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are  printed.   The  definition  of  the  dimens ionless  amplitude  and 
phase  angles  are  given  on  the  print-out  and  will  be  given  here 
also  for  completeness: 

PI  =  pressure  (half)  amplitude  in  height  of  water  column/ (ha  1 f ) 

amplitude  of  the  motion  in  surge. 
P2  =  pressure  (half)  amplitude  in  height  of  water  column/ (hal f ) 

amplitude  of  the  motion  in  heave. 
P3  =  pressure  (half)  amplitude  in  height  of  water  column/ 

(half)  amplitude  of  the  motion  in  sway. 
P4  =  pressure  (half)  amplitude  in  height  of  water  column/ 

(half)  amplitude  of  the  motion  in  roll  times  ABAR. 
P5  =  pressure  (half)  amplitude  in  height  of  water  column/ 

(half)  amplitude  of  the  motion  in  yaw  times  ABAR. 
P6  =  pressure  (half)  amplitude  in  height  of  water  column/ 
(half)  amplitude  of  the  motion  in  pitch  times  ABAR. 
These  definitions  are  the  same  as  stated  in  Eqs.  (2.20). 

The  actual  pressures  at  any  given  node  point  could  be 
expr es sed  as  : 

P   =  PI  p gX  °  a  cos  (PHI  -  at) 
?2    -  P2  p  gX2°  a  cos  (PH2  -  at) 


o  - 


P   =  P3  p  gX    a  cos  (PH3  -  at) 


o  - 


P,  =  P4  p  gX    a  cos  (PH4  -  at) 


o  - 


P   =  P5  p  gX    a  cos  (PH5  -  at) 


o  - 


P,  =  P6  p  gX    a  cos  (PH6  -  at) 
Py  =  P7  p  g  (H/2)  cos  (PH7  -  at) 
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where  X    denote  the  dimens ionles s  half  amplitudes  of  the  motion 
i 

as  defined  in  Eq .  (2.1)  where 

X^  =  X^/a 

o 
X°2  =  X2  /a 

x°3  =  x3°/i 

Since  the  motion  of  the  body  is  defined  by 

X. (t)  =  X°.  cos ( 6  -at) 
1         l       i 

it  is  clear  the  positive  phase  angles  represent  a  lag  of  the 
pressure  with  respect  to  the  displacement  in  any  given  degree  of 
freedom . 

The  next  block  of  output  represents  the  three  components  of 
fluid  velocity  (in  body  axes)  associated  with  motion  in  the  six 
degrees  of  freedom  (1-6),  and  wave  interaction  with  the  fixed  hull, 
(7  ), cor responding  to  each  node  on  the  complete  immersed  surface. 
(Note:   if  NBA  =  1  this  information  is  not  computed  and  the 
output  is  suppressed). 

The  velocity  component  dimens ionless  amplitudes  are  defined, 
for  example ,  as 

VX  =  velocity  (half)  amplitude  in  the  x ' -direction/ 
velocity  (half)  amplitude  of  the  body. 

The  phase  angles  are  defined  with  respect  to  the  displacement 
of  the  body  (positive  lag).   The  x'-component  of  velocity  could 
be  expressed,  for  example,  as 
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surge  (1)  V  =VX  X    ao    cos(pxl-  at) 

X  J. 

heave  (2)  V  =VX  X°   ao  cos(PX2-  at) 

x  z. 


e tc  . 

The  MASS  MATRIX  is  next  printed  exactly  as  read-in  through 

2  2 

data  cards.   The  units  are  slugs  and  slugs-ft  ,  or  kg  and  kg-m  . 

The  SPRING  CONSTANT  MATRIX  is  printed  next  also  as  it  was 
read  in  through  the  six  data  cards.   Note  that  for  a  free-floating 
body  all  of  the  values  in  this  array  are  zero. 

The  next  set  of  data  represents  the  dimens ionless  wave  load 
coefficients.   The  x,y,z  subscripts  represent  body  axes  so  FX,  for 
example,  represents  the  dimens ionless  excitation  force  in  the 
x'  direction.   The  definition  of  the  coefficients  and  phase 
angles  are  given  on  the  print-out.   The  actual  wave  force  and 
moments  can  be  expressed  as  a  function  of  time  as: 

F   =  FX  pga3  (H/2a)cos  (PHASE  FX-ot) 

X 

F       =    FY    pga3(H/2a)     cos     (PHASE    FY-ot) 

y 

F   =  FZ  pga3(H/2a)  cos  (PHASE  FZ-at) 

M   =  MX  p ga  (H/2a)  cos  (PHASE  MX-ot) 
x 

M   =  MY  p ga  (H/2a)  cos  (PHASE  MY-at) 

4  - 
M   -  MZ  p ga  (H/2a)  cos  (PHASE  MZ-ot) 

Note:   the  wave  crest  passes  the  origin  of  the  coordinate  system 
at  t  =  0.   Thus,  a  positive  phase  angle  represents  a  lag  with 
respect  to  the  time  the  wave  crest  passes  the  origin. 

The  HYDRODYNAMIC  ADDED  MASS  COEFFICIENT  MATRIX  is  next 
printed  followed  by  the  HYDRODYNAMIC  DAMPING  COEFFICIENT  MATRIX. 
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These  parameters  are  defined  through  Eqs.  (4.9  -  4.13).   However, 
to  emphasize  their  relationship  to  the  velocity  and  acceleration 
their  definitions  will  be  given  here  again. 

The  force  (or  moment)  in  the  i-direction  (in  body  axes)  due 
to  motion  in  the  i-direction  is  defined  as  F. .  and  is  given  in 

terms  of  the  added  mass  and  damping  coefficients,  M. .  and  N. . ,  as 

_4     ••        _A 
F..  =  -(1.0  or  a)  {p a  M..  X.  +Doa   N..  X.} 

where  1.0  is  used  when  F. .  represents  a  force  (i  =  1,2,3)  and 

a  is  used  when  F. .  denotes  a  moment  (i  =  4,5,6) .   In  the  above 

expression 

X.  =  acceleration  in  i-direction  (Note,  X.  =  X./a 
J  3     J 

for  j  =  1,2,3  and  X.  =  0.  for  j  =  4,5,6). 

X.  =  velocity  in  i-direction 
3 

The  next  block  of  output  is  labeled  HASKINDS  RELATIONS. 
Here  the  excitation  forces,   moments,  phase  angles  are  computed 
by  use  of: 

(1)  the  near  field  pressure  distribution 

(2)  the  far  field  (radiation)  potential  by  use  of  Haskind's 
relations 

and  the  damping  coefficients  are  computed  by 

(1)  the  near  field  pressure  distribution 

(2)  the  energy  radiated  at  infinity. 

The  results  computed  by  (1)  and  (2)  are  then  compared  and  a 
percent  difference  computed. 

The  use  of  the  Haskind's  relations  and  energy  balance 

represents  a  self-check  on  the  accuracy  of  the  results  computed 
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by  the  pressure  distribution.   The  results  computed  by  use  of  the 
pressure  distribution  should  be  more  accurate  than  those  computed 
from  the  far  field  behavior  of  the  potential  but  the  two  results 
should  be  fairly  close.   Large  percentage  differences  indicate 
that  there  may  be  some  problem  and  the  solution  may  not  be  valid. 
However,  if  the  values  of  the  parameters  being  compared  are  very 
smalljlarge  percentage  differences  should  be  expected  and  should 
cause  no  significant  error  in  the  final  results. 

Finally  the  DYNAMIC  RESPONSE  results  are  printed.   The 
response  amplitude  and  phase  angle  in  the  six  degrees  of  freedom 
are  printed.   The  amplitudes  are  defined  as: 

Surge:   Half  amplitude  of  the  surge  motion/wave  (half)  amplitude 
Heave:   Half  amplitude  of  the  heave  motion/wave  (half)  amplitude 
Sway:   Half  amplitude  of  the  sway  motion/wave  (half)  amplitude 
Roll:   Half  amplitude  of  the  roll  motion  in  radians  x  ABAR/wave 

(half)  amplitude. 
Yaw:    Half  amplitude  of  the  yaw  motion  in  radians  x  ABAR/wave 

(half)  amplitude. 
Pitch:   Half  amplitude  of  the  pitch  motion  in  radians  x  ABAR/wave 
(half)  amplitude. 

The  phase  angles  of  the  dynamic  response  are  defined  with 
respect  to  the  crest  of  the  incident  wave  in  radians.   The  crest 
of  the  incident  wave  passes  the  origin  at  t=0  so  a  positive  phase 
angle  indicates  that  the  displacement  of  the  body  lags  the  crest. 
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9.   EXAMPLE  COMPUTATIONS 

In  this  chapter  a  simple  example  calculation  is  carried  out 
for  a  rectuangular  barge  (box)  which  is  90  x  90  meters  in  plan 
with  a  draft  of  40  meters.   Results  are  computed  for  a  series  of 
different  wave  periods. 

Hull  Corner  Node  Points 

In  Figure  9  is  shown  the  subdivision  of  the  hull;  since  the 
immersed  surface  possesses  double  symmetry  it  is  necessary  to 
describe  only  one  quarter  of  it  by  use  of  data  cards. 

The  subdivision  indicated  in  Figure  9  is  rather  coarse. 
There  are  12  panels  on  the  quarter-body  and  48  panels  all  together 
on  the  total  immersed  surface.   The  quarter-body  is  described  by 
use  of  the  19  corner  node  points  so, 

NC  =  19 

NA  =  12 

NB  =  48 

for  the  immersed  surface  described  in  Figure  9. 

The  immersed  surface  was  input  with  an  input  coordinate 

sys  tern  placed  relative  to  the  body  as  indicated  in  Figure  9 

with  the  bottom  represented  by  the  y"  =  0  plane.   Since  the  draft 

was  40.0  m  the  inertial  coordinate  system  was  placed  with  origin 

at  the  position  x"  =  z"  =  0  and  y"  =  40.0.   This  is  input  to  the 

computer  program  by  use  of  a  data  statement  specifying  the 

parameters  as  follows: 

XREF  =  0.0 
YREF  =  40.0 
ZREF  =  0.0 
ANGREF  =  0.0 
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The  computer  "echo  check"  on  the  input  data  is  shown  in 
Table  2.   The  location  of  the  inertial  reference  frame  in  the  input 
reference  frame  is  first  printed  in  the  table.  The  rotation  angle 
is  zero . 

TA3LE  2 


DOTATION  ANGLf:  = 
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The  coordinates  of  the  corners  of  the  panels  are  printed 
under  the  heading:   COORDINATE  OF  PANEL  CORNERS.   A  comparison 
of  the  numbers  listed  here  can  be  seen  to  agree  with  the  location 
of  the  corners  (in  the  inertial  reference  frame)  shown  in  Figure  9. 

Next,  under  the  heading  CORRESPONDENCE  TABLE  can  be  found 
the  information  which  describes  the  given  panels  (numbered  1-12) 
in  terms  of  the  corner  node  point  indices.   Under  NODE  1,  NODE  2, 
etc.  are  listed  the  four  corner  node  points  which  go  to  make  up 
a  given  panel  listed  in  clockwise  order  when  observing  from  the 
"water  side"  of  the  panel.   A  comparison  of  these  results  with 
Figure  9  will  confirm  this. 

Table  3  gives  the  computer  print-out  of  the  location  of  the 
panel  node  points,  the  components  of  the  outward  unit  normal 
vectors  at  the  immersed  surface  (in  body  axes)  and  the  areas  of 
the  panels.   The  node  points  are  calculated  by  the  computer,  as 
described  in  Appendix  A,  and  represent  the  centroids  of  the 
panels . 

In  Table  4  is  printed  the  volume  and  surface  area  of  the 

immersed  surface  as  computed  on  the  basis  of  the  information  given 

in  Table  3.   The  volume  is  computed  by  three  different  methods, 

i.e.,   by 

n 

(x-proj)   i^1   i  x±  l 


(y-proj)  =   X!   y-  n    AS 


i-1   X   yi 
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V,       .x  =   Z]   z  .  n    AS  . 
(z-proi )     ,  -   1   z  .    1 

The  volume  calculations  in  this  way  use  all  of  the  information 
in  Table  3  and  if  the  volumes  computed  by  the  different  methods 
are  equal  to,  say,  5  or  6  decimal  places,  it  is  quite  likely 
that  the  input  data  is  correct.   This  represents  a  simple  quick 
verification  of  the  input  data  but  cannot  replace  the  Calcomp 
plot  for  ease  and  reliability. 

The  different  indices  assigned  to  the  dif  f er ent  levels  of 
the  panel  node  points  are  also  given  in  Table  5.   It  may  be  seen 
by  looking  over  the  results  printed  in  Table  3  that  the  node 
points  fall  at  only  three  different  levels:   -10.0,  -30.0  and 
40.0.   As  indicated  in  Table  5,  and  as  can  be  verified  by  Figure 
9,  panel  nodes  1-4  are  on  the  bottom  and  are  located  at  the  -40.0 
level  (LEVEL  NO  1),  panel  nodes  5,6,9,10  are  located  at  the  -10.0 
level  (LEVEL  NO  2),  and  panel  nodes  7,8,11,12  are  at  the  -30.0 
level  (LEVEL  NO  3) . 

The  next  block  of  computer  print-out  gives  the  WAVE  NUMBER 
and  lists  the  parameters,  a  =  2iTa/L,  v  =  a  tanh(ah)  ,  where 
h  =  h/a,  water  depth,  wave  period,  angle  of  incidence  and  the  value 
of  NB. 

Following  this  information,  the  computer  print-out  contains 
the  pressure  distribution  and  velocity  distribution  as  previously 
described.   After  this  an  "echo  check"  on  the  MASS  MATRIX  and 
SPRING  CONSTANT  MATRIX  is  printed  in  the  same  units  as  input  by 
way  of  the  data  cards.   This  information  appears  to  be  straight- 
forward and  will  not  be  given  here. 
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The  next  block  of  information  on  the  print-out  is  the  dimen- 
sionless  excitation  force  and  moment  coefficients.   For  the  sample 
run  this  is  presented  in  Table  6.   It  may  be  noted  that  the  angle 
of  incidence  of  the  wave  was  zero  so  naturally  Table  6  indicates 
FZ  =  MX  =  MY  =  0.0. 

Table  7  contains  the  next  three  blocks  of  computer  print- 
out.  These  three  blocks  represent  the  dynamic  and  hydrostatic 
force  coefficient  associated  with  the  motion  of  the  body.   In 
making  this  example  calculation,  NSS  was  set  to  1  even  though  it 
could  have  been  set  to  2  since  the  body  is  symmetrical  about  the 
x'-y'  plane  and  the  angle  of  incidence  was  zero.   Since  NSS  =  1 
all  six  degrees  of  freedom  of  the  body  were  considered  and,  accor- 
dingly, all  of  the  added  mass  and  damping  coefficients  were  com- 
puted.  Had  NSS  been  set  to  2  only  (M.,,  M.„  and  M  .,)  and 

ll    i2        16 

(N.,,  N.„  and  N.,)  would  have  been  computed.   The  value  of  the 

ll     i2         10 

other  coefficients  would  be  meaningless. 

It  may  be  noted  that  the  surge  and  sway  added  mass  and 
damping  coefficients,  M    and  M„„,  and  N  -  and  N~„,  are  equal. 
This  is,  of  course,  due  to  the  fact  that  the  body  is  symmetrical 
in  plan  view. 

It  is  also  interesting  to  note  that  the  roll,  yaw  and  pitch 

damping  moments,  i.e.,  N   ,  N    and  N,,,  are  quite  small.   This, 

4  4     5  5        do 

in  fact,  appears  to  be  typical  of  most  cases,  and,  therefore, 
must  be  taken  into  consideration  when  interpreting  that  data. 
It  may  be  recalled  that  the  computed  damping  is  wave  damping 
only  and  this  tends  to  be  small  in  the  angular  modes.   In  such 
cases  viscous  damping  may  become  important  and  it  may  be  worthwhile 
to  estimate  the  viscous  damping. 
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TABLE   6 
****  yAVE  LaA0  COEFFlCItNTS  AND  PHASE  ANGLES  FOR  THE  CAISSON  **** 

FX,  FY,  FZ     =  AMPLITUDE  GF  FC*C '  /  >•>  FC*G*  ABA  R**3  (  H/2  +  AEAR  ) 

NX,  MY,  MZ    =  AHPLITUOE  GF  HQnENT  /RhQ  *G  *ABAR  **4  (  H/  2*  Ab  AR  I 

PHASE  «NGLES  =  IN  RADIANS  MEASURE 0  wITF  RESPECT  TO  TlPf- 
WAVl  CREST  IS  A^  JUGIN  (POSITIVE  =  LAG) 

PHASE  FX=  -  1.  3296 

PHASE  FY=  -0.<i^Z<* 

PHASE  FZ^   0.5Cld 

PHASE  MX=   3.1C75 

PHASE  MY  =   2.2377 

PHASE  MZ  =  -1.2726 

TABLE   7 

****  HYDRQQYNAMIC    ADDED    MASS  CnrFFICIcM  MATRIX    **** 

2.7645  -J.0000  -0.0000  -0.0300  -C.CGCO  C.2076 

C.OOOO  2.3531  -0,0000  0.0000  0.0000  -C.0000 

C.OOOO  0.0000  2./6h5  -0.2076  C.OOGO  C.OOOO 

C.OOOO  J. 0000  -0.2223  J. 4534  0.0000  C.OOOO 

C.OOOO  0.0000                  0.0000  -0.0000  C.6559  -C.OOOO 

0.2224  -0.0000                  J. 0000  0.0000  -O.OOCO  0.4534 


FX  = 

2.9564 

FY  = 

1.0021 

FZ  = 

0.0000 

MX  = 

0.0  000 

MY  = 

0.0000 
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The  next  block  of  information  printed  in  Table  8  is  the 
Energy  and  Haskind's  relations  check  on  the  accuracy  of  the 
solution.   Here  the  excitation  forces,  moments  and  phase  angles 
as  well  as  the  diagonal  of  the  damping  matrix  is  printed.   The 
calculation  of  these  parameters  is  carried  out  both  through  a 
pressure  integration  over  the  immersed  surface  and  using  the 
"far-field"  potential  through  Haskind's  relations  and  the  energy 
balance.   The  two  sets  of  values  and  the  percent  difference 
between  the  two  results  are  printed. 

It  should  be  kept  in  mind  that  in  this  run  FZ  =  MX  =  MY  =  0 
so  the  percent  differences  in  these  cases  is  meaningless.   Gen- 
erally, even  for  this  very  coarse  grid  the  two  sets  of  results 
agree  rather  well,  at  least  when  the  values  are  not  extremely 
small.   For  instance,  FX  and  FY  differ  by  only  5  to  6  percent. 
MZ,  however,  is  rather  small  and  the  two  results  differ  by  about 
33  percent.   This  is,  generally,  the  case;  the  moment  tends  to 
be  rather  sensitive  to  grid  size  and  it  requires  a  finer  grid 
to  get  accurate  moment  results. 

The  values  of  N1 1 ,  N99  and  N«»  differ  by  less  than  6  percent 
when  computed  by  the  two  different  methods.   However,  the  damping 
in  pitch,  yaw  and  roll  are  very  small  and  tend  to  differ  by  a 
sizable  amount.   The  difference  is,  however,  not  too  important 
because  such  small  values  of  the  damping  coefficient  would  not 
affect  the  response  results  except  at  resonance. 
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The  dynamic  response  of  the  floating  body  is  given  in  Table 
9.   These  results  indicate  that  the  surge  (half -amplitude )  is 
•5565(H/2)  where  H  denotes  the  wave  height  (trough  to  crest).   In 
this  particular  run  the  heave  response  shows  a  little  dynamic 
amplification  since  the  dimens ionless  response  is  slightly  greater 
than  1.0.   The  heave  half -amplitude  is  1.2(H/2).   The  pitch  half- 
amplitude  is  .0758  x  H/(2a),  in  radians,  where  a  =  45.0  m  in  this 
example  problem. 


TABLE  9 
DYNAMIC  RESPONSE 

Surge  (X)  **  X/ETA  =     0.55650  PHASE  =   1.51353 

Heave  (Y)  **  Y/ETA  =     1.19997  PHASE  =   2.22306 

Sway   (Z)  **  Z/ETA  =     0.00000  PHASE  =  -1.11987 

Roll   (X)  **  X/ETA  =     0.00001  PHASE  =  -0.71514 

Yaw    (Y)  **  Y/ETA  =     0.00000  PHASE  =   0.17472 

Pitch  (Y)  **  Y/ETA  =     0.04758  PHASE  =  -1.61674 


102 


10.   NUMERICAL  RESULTS 

The  purpose  of  this  chapter  is  to  present  some  typical 
numerical  results  for  several  simple  geometric  configurations 
and  to  make  comparisons  with  experimental  results.   It  appears 
that  very  little  experimental  data  is  in  fact  available  for  com- 
parison with  the  computed  results  but  in  cases  where  data  is 
available  comparisons  are  made.   In  addition,  a  comparison  of 
the  effect  of  subdivision  size  on  typical  results  are  presented. 
Numerical  Results  for  a  90m  x  90m  x  40  m  Free-Float ing  Box: 

Computer  calculations  have  been  carried  out  for  the  case 
of  the  floating  box,  90  meters  by  90  meters  in  plan  by  40  meters 
draft,  which  was  discussed  in  Chapter  10.   The  hydrodynamic  para- 
meters as  well  as  the  dynamic  response  to  wave  motion  are  plotted 
as  a  function  of  wave  period  for  two  different  grid  sizes, 
NB  =  48  and  108. 

The  two  different  subdivisions  of  the  floating  box  are 
shown  in  Figures  10  and  11  for  the  NB  =  48  and  108  panel  layouts. 
It  may  be  noted  that  the  figures  show  only  one  quarter  of  the  box 
because  the  box  has  two  planes  of  symmetry  and,  accordingly,  in 
the  two  cases  NA  =  12  and  27. 

Model  tests  were  carried  oat  for  the  configuration  in  ques- 
tion at  the  Vassdrags  Og  Havnelaboratoriet  (Rivers  and  Harbors 
Laboratory)  in  Trondheim,  Norway.   These  results  which  were  re- 
ported by  Faltinsen  and  Michelsen  (11)  are  compared  with  the 
numerical  results  presented  here. 
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40m 


FIGURE   10:    FLOATING  BOX,    90M  x  90M  x  40M.    NB=48  PANELS 
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FIGURE  11:    FLOATING  BOX,    90M  x  90M  x  40M.    NB=108  PANELS 
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The  surge  and  heave  excitation  forces  are  shown  in  Figures 
12  and  13,  respectively.   The  fact  that  there  is  very  little 
difference  between  the  two  sets  of  results  corresponding  to 
NB  =  48  and  108  indicates  that  the  solution  tends  to  converge 
very  rapidly  in  increasing  the  number  of  panels  (or  decreasing 
the  grid  size),  at  least  for  the  case  of  the  excitation  forces. 
It  might  be  noted  that  slower  convergence  should  be  expected  at 
small  periods . 

The  pitch  excitation  moment  was  computed  also.   However, 
as  it  turned  out  in  this  particular  case,  the  moment  with  respect 
to  the  body  axes  (i.e.,  the  centroidal  axes  was  very  nearly  zero 
throughout  the  complete  frequency  range  and  the  results  were, 
therefore,  masked  by  the  "noise"  caused  by  the  numerical  error.) 
For  all  practical  purposes  the  pitch  excitation  moment  was  zero. 

The  surge  added  mass  coefficient  is  shown  in  Figure  14,  and 
here  again  the  results  corresponding  to  NB  =  48  and  108  indicate 
that  the  solution  has  converged.   The  experimental  results  pre- 
sented for  comparison  appear  to  agree  rather  well  with  the  theory 

The  corresponding  damping  coefficient  in  surge  is  shown  in 
Figure  15.   These  results  are  typical  of  most  damping  results. 
At  high-frequency  and  low- frequency  oscillation,  the  wave-making 
damping  tends  to  vanish  since  no  waves  are  generated  at  the  two 
extremes.   In  the  intermediate  frequency  range,  however,  the 
wave-making  reaches  a  maximum  and  the  damping  coefficient  like- 
wise shows  a  maximum.   The  experimental  results  shown  on  the 
figure  appear  to  agree  fairly  well  with  the  theory. 
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The  added  mass  coefficient  in  heave,  M   >  is  shown  in  Figure 
16  along  with  experimental  values.   These  results  show  very  little 
frequency  dependence,  particularly  when  compared  to  the  surge 
damping  coefficient,  N   .   This  results  from  the  fact  that  the 
surface  primarily  involved  in  heave  force  is  the  bottom  and  this 
surface  is  fairly  distant  from  the  free  surface.   Thus,  free 
surface  effects  which  are  frequency  dependent  are  not  pronounced 
in  the  case  of  heave  motion  whereas  they  are  much  more  pronounced 
in  the  case  of  surge  where  surfaces  near  the  free  surface  are 
involved . 

The  computed  heave  damping  coefficient  is  shown  in  Figure  17 
along  with  experimental  values.  It  is  interesting  to  note  that 
compared  to  the  surge  damping  coefficient  the  heave  damping  co- 
efficient is  rather  small.  This  results  from  the  fact  that  the 
bottom  surface  which  acts  as  the  primary  wave-maker  in  heave  is 
rather  well  removed  from  the  free  surface  and  is,  therefore,  not 
too  effective  in  generating  waves  and  resulting  damping. 

It  is  also  of  interest  to  note  that  the  experimental  values 
of  N    fall  well  above  the  theoretical  results.   This  difference 
is  attributed  to  viscous  damping  and  appears  to  be  pronounced 
in  a  relative  sense  because  the  computed  wave-making  damping  is 
small . 

The  added  moment  of  inertia  coefficient  shown  in  Figure  18 
is  essentially  independent  of  frequency  (period).   This  appears 
to  be  typical  of  all  results  associated  with  caissons  of  the  type 
in  question.   Moreover,  the  damping  coefficient  in  pitch,  N^, 
which  is  not  plotted,  tends  to  be  extremely  small.   These  two 
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features  result  from  the  fact  that  frequency  dependence  and 
damping  is  a  result  of  wave  production  and  very  little  wave  gen- 
eration occurs  as  a  result  of  pitching  motion.   This  is  due  to  the 
fact  that  the  bottom  surface  is  far  removed  from  the  free  surface 
and  the  vertical  sides  do  not  move  in  the  normal  direction  very 
much  as  a  result  of  pitching  motion. 

The  added  mass  coefficient  in  yaw,  M  ,  is  shown  in  Figure 
19  along  with  experimental  results.  The  agreement  here  appears 
to  be  generally  good. 

Finally,  the  dynamic  response  of  the  box  free  floating  in 
waves  is  shown  in  Figure  20.   The  pitch  response , 0, / (H/ 2a ) ,  is 

D 

extremely  small  as  was  expected.   The  natural  frequency  in  heave 
occurs  at  about  16.5  seconds  period  and,  therefore,  considerable 
dynamic  amplification  occurs.   In  the  case  of  the  heave  response 
several  experimental  values  are  shown  which  tend  to  agree  well 
with  the  theory. 
Shallow-Draft  Barge 

The  dynamic  heave  and  pitch  response  has  also  been  computed 
for  comparison  with  experimental  results  for  a  shallow- draf t 
barge.   One  quarter  of  the  immersed  surface  of  the  barge  is  shown 
in  Figure  21  as  it  was  subdivided  for  numerical  computations.   It 
may  be  noted  that  since  the  draft  was  rather  small  the  source 
strength  might  be  expected  to  show  considerable  variation  in 
the  vertical  direction  and  around  the  corner.   However,  in  the 
horizontal  direction  the  source  strength  should  vary  rather  slowly 
Accordingly,  the  panels  were  constructed  with  a  rather  large 
aspect  ratio  as  indicated  in  Figure  21. 
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Calculations  have  been  made  to  compare  with  the  wave  channel 
test  results  of  Kim,  Henry  and  Chou  (17)  for  a  shallow  draft  barge 
of  1.704  feet  beam,  by  2.558  feet  length,  by  0.16  feet  draft. 
The  location  of  the  center  of  gravity  and  pitch  gyradius  is  shown 
in  Figure  22. 

The  theoretical  heave  and  pitch  response  is  compared  with 
the  experimental  values  in  Figure  22.   The  results  show  excellent 
agreement  between  the  theory  and  experiment  for  the  case  of  heave 
motion.   This  apparently  results  from  the  fact  that  the  wave -making 
damping  in  heave  for  a  shallow-draft  barge  is  substantial  and, 
therefore,  the  heave  results  are  not  too  sensitive  to  viscous 
damping . 

The  pitch  results  also  show  fairly  good  agreement.    While 
there  appears  to  be  some  scatter  in  the  pitch  results,  it  appears 
that  near  resonance  the  experimental  values  of  the  response  fall 
below  the  theory.   Here  again  this  is  attributed  to  viscous 
effects  which  tend  to  be  significant  compared  to  the  rather  small 
amount  of  wave-making  damping. 
Disc  Buoy 

As  a  further  example,  the  dynamic  heave  and  pitch  response 
for  a  rather  shallow-draft  disc-buoy  is  shown  in  Figure  23.   Un- 
like some  of  the  previous  results,  the  pitch  response  as  well  as 
the  heave  response  shows  excellent  agreement  between  the  theory 
and  experiment.   This  is  likely  due  to  the  fact  that  the  corners 
on  the  disc  buoy  configuration  are  beveled  rather  than  sharp.   This 
probably  reduces  some  of  the  separation  and  resulting  viscous 
damping.   Consequently,  the  experimental  results  are  in  good 
agreement  with  the  calculations. 
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Floating  Semi-submerged  Sphere 

As  an  example  showing  the  effect  of  finite  depth  on  a  floating 
body,  numerical  results  are  presented  for  a  semi-submerged  sphere 
floating  in  waves.   Figures  24,  25  and  26  show  the  excitation 
forces,  added  mass  and  damping  in  heave,  and  added  mass  and 
damping  in  surge,  respectively.   In  general,  the  effect  of  the 
finite  depth  tends  to  be  more  pronounced  at  low- frequency  than  at 
high-frequency. 

Using  the  results  of  Figures  24-26  the  equations  of  motion 
for  the  free  floating  sphere  give  the  heave  and  surge  response 
presented  in  Figures  27  and  28.   The  amplitude  ratios  are  shown 
in  Figure  27  and  the  phase  shift  angles  of  the  response  measured 
with  respect  to  the  phase  of  the  wave  crest  at  the  center  of  the 
body  is  shown  in  Figure  28.   In  general,  there  appears  to  be  little 
effect  on  the  finite  depth  for  the  sphere  configuration. 
Floating  Vertical  Cylinder 

The  added  mass  and  damping  coefficients  for  the  circular 
caisson  configuration  are  shown  in  Figure  29  as  a  function  of 
the  frequency  parameter  a   a/g  wherein  the  depth  ratio  h  =  h/a, 
represents  a  parameter.   The  horizontal  force  coefficient  is  shown 
in  Figure  30.   These  results  for  the  horizontal  force  coefficient 
agree  well  with  the  theoretical  results  of  Garrett  (18). 

The  dynamic  response  in  surge,  heave  and  pitch  is  presented 
in  Figures  31,  32  and  33,  respectively.   In  general,  the  results 
appear  typical  of  deep  water  results  and  little  effect  of  water 
depth  is  evident.   The  pitch  response  shown  in  Figure  33  shows 
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the  large  dynamic  amplification  at  resonance  which  tends  to  be 
typical  of  essentially  all  floating  caissons  and,  here  again 
results  from  the  rather  small  wave-making  damping.   This  high 
resonant  peak  should  not  be  expected  in  practice,  however,  due 
to  the  fact  that  viscous  effects  are  undoubtedly  considerably 
larger  than  the  pitch  wave-making  damping  used  to  compute  the 
results  shown  in  Figure  33. 
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CjJ=  FL(max)/pga  (H/2) 
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FIGURE  24   WAVE  EXCITATION  FORCES  ON  A  SEMI-SUBMERGED  SPHERE 
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FIGURE  25   ADDED  MASS  AND  DAMPING  IN  HEAVE  FOR  A  SEMI 
SUBMERGED  SPHERE. 
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FIGURE  27  HEAVE  AND  SURGE  RESPONSE  OF  A  SEHI-SUBI1ERGED  SPHERE 
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HEAVE,  \\f 


FIGURE  28  PHASE  SHIFT  ANGLES  OF  THE  HEAVE  AND  SURGE  RESPONSE  OF 
SEIII-SUBME^GER  SPHERE. 
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FIGU-JE  29  ADDED  MASS  AND  DAMPING  COEFFICIENTS  FOR  A  FLOATING 
VERTICAL  CIRCULAR  CYLINDER. 
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FIGURE  30  HORIZONTAL  FORCE  COEFFICIENT  FOR  A  FLOATING  VERTICAL 
CIRCULAR  CYLINDER. 
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FIGURE  31  SURGE  RESPONSE  OF  A  FLOATING  VERTICAL  CIRCULAR  CYLINDER 
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FIGURE  32  HEAVE  RESPONSE  OF  A  FLOATING  VERTICAL  CIRCULAR  CYLINDER. 
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FIGURE  33    PITCH  RESPONSE  OF  A  FLOATING  VERTICAL  CIRCULAR  CYLINDER 
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APPENDIX  A: 

DETERMINATION  OF  THE  CENTROIDAL  LOCATION,  AREA  AND  UNIT  NORMAL 

VECTOR  FOR  A  PANEL 

The  computer  program  is  designed  to  accept  the  coordinates 
of  the  corners  of  the  facets  as  inputs  and,  in  addition,  a  set 
of  cards  indicating  which  corners  make-up  a  given  panel.   A 
panel  is  described  by  the  four  corner  indexes  read 
clockwise  moving  around  the  panel  in  such  a  way  that  the  panel 
is  always  on  the  right.   For  example,  panel  5,  as  indicated  in 
Figure  1A,  is  specified  by  a  data  card 


FIGURE  1A 


of  the  following  form: 

Panel    No.  Corner    //l  Corner   ill        Corner    //3  Corner    //4 


25 


31 


50 


18 


In  working  with  an  individual  panel  the  first  corner  (25)  in 
the  above  example   is  designated  as  1,  the  second  corner,  moving 
around  the  panel  keeping  the  inside  of  the  panel  on  the  right, 
is  labled  2,  etc.   The  panel  is  re-labled  as  shown  in  Figure  2A 
with  indices  1-4. 

The  area  of  the  panel  is  computed  by  first  dividing  it  into 
two  triangles  by  connecting  corner  2  with  the  diagonal  corner  4. 
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(<*.A,,2,) 


FIGURE  2A 


The  area  of  triangle  1-2-4  is  then  computed  by  use  of  a  standard 
formula , 


where 


s  =  i  (dlt  +  d„.  ♦  a4l) 


and 


d„  = 

dz.4.  = 
<L.   = 


vat- 

-*r  + 

^y«- 

y,;z 

-h 

(**- 

-z,;z 

V(x4- 

-  x)1  + 

(* 

-y.f 

+- 

(*»- 

■zo*- 

V(^~  x.f  +  Cy.-yJ2  i-  Cz,-  z4f 


The  centroid  of  a  triangle  lies  at  the  average  of  the  co- 
ordinates describing  its  corners  as, 


3 


Z, 


>2-4 


=   4r(Z.  t  2-z.  +Z-^) 
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The  area  and  centroid  location  of  triangle  234  is  computed 
in  a  similar  manner.   The  total  area  for  the  panel  is  then  simply 


A  -   A 


12.4- 


Ar 


34- 


and  coordinates  of  the  centroid  of  the  panel  is  obtained  from 


It  is  also  necessary  to  determine  the  x , y , z-components  of 
the  unit  vector  normal  to  the  panel  and  directed  outward  into 
the  fluid.   This  is  accomplished  by  defining  the  vectors,  T} 
and  *7J   ,  which  extend  from  point  1  to  3  and  from  point  2  to  4  , 
respectively,  as  indicated  in  Figure  3A.   The  corners  are  numbered 


FIGURE  3A 

clockwise  when  viewed  from  the  outside  (water  side)  of  the  caisson 
and,  therefore,  according  to  the  right  hand  rule  the  cross  product 

represents  a  vector  having  direction  perpendicular  to  the  plane 
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1-2-3-4  and  pointing  outward  into  the  fluid.   The  unit 
normal  vector  is,  then,  obtained  by  normalyzing  this  vector  to 
unit  length, 


w 


here  |ftl  =  |  and  R  ^  Lfl^  jH^j  fe  fl^   ,  the  components   Hx , 
flu  and  rifc  represent  the  direction  cosines. 
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APPENDIX  B: 

INTEGRATION  OF  1/R  AND  DERIVATIVES  OF  1/R  OVER  A  PANEL 

In  order  to  evaluate  the  elements  of  the  matrices  Q(   and  (3 
by  use  of  the  definitions  given  in  equations  (  5. 3~)  and  (5.7) 
the  necessity  arises  to  evaluate  integrals  of  1/R  and  a/^bx, 

0/dL\    and  o/3^  of  1/R.   In  this  appendix  the  method  of 
integration  is  outlined. 

In  order  to  carry  out  the  integrations,  a  local  coordinate 
system  is  defined  in  such  a  way  that  the  panel  lies  in  the  x-y 
plane  of  the  local  coordinate  system  and  the  unit  normal  vector 
is  parallel  to  the  z  axis.   For  definiteness  the  coordinate 
system  is  attached  to  the  element  in  such  a  way  that  point  1 
is  located  at  the  origin  and  point  4  lies  on  the  x-axis. 
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FIGURE  IB 


Consider  first  the  integral, 


<t> 


-kds 


cU_jL 


AS 


AS 


(IB) 
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Following  a  procedure  similar  to  Hess  and  Smith (16)  for  the 
velocity  components,  Faltinsen  and  Michelsen(  11 )  integrated 
equation  (IB   with  respect  to  ^  to  obtain 


h 


(2B) 


where 


(3B) 


The  integrals  in  (2B)  may  be  evaluated  through  numerical  integ- 
ration in  most  cases. 

However,  it  may  be  noted  that  the  integrand  of  the  first 
integral  is  singular  when  z  =  O,     5=  x  and  y  -  ?   <  0.   In 


r*» 


this  case  the  integral  may  be  replaced  by 
St. 

A 


(4B) 


-fj-f-u-ij  *■  WiJ^1^71^} 


ds 
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The  first  integral  in  (4B)  can  be  integrated  analytically,  the 
result  of  which  is 

(5B) 


The  second  integral  in  (4B)  has  no  singularity  in  the  integrand 
and  consequently  presents  no  dif f iculat ies  in  numerical  integ- 
ration.  Difficulties  with  singularities  in  the  integrand  of 
the  other  integrals  in  equation  (2B)  are  handled  in  a  similar 
manner . 

In  the  evaluation  of  CV^;  as  defined  by  equation  (4.5)  or 
yC..        .      Yu .  ■        and  X  .  .  as  given  by  equations  (4.11)  it  is 
necessary  to  evaluate  integrals  of  the  form: 


AS 


(6B) 


(7B) 


*-SL 


c/S 


AS  «*  (8B) 


where 

R=  V  (X-sr  i-   Cy-j^  +   **-' 

Equations  (6B-8B)  have  been  integrated  by  Hess  and  Smith3 
the  results  of  which  are  summarized  as: 
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OB) 


-h 


d 


3+ 


[*3    +A-4.     +d^    )  d+,  K/^+h,**** 


5,-  5 


(10B) 


+  J,-  -^V  in  /*j  *  <*+-d3+ )  +  g+-X,Jn  A^A-  *?+,   ) 


-      /a 


r2 


4>  =   fan  I  ^^  e,-  h,\ 


+   hn    fmu£jLJL±±)  _     Ian"  / 


2  rLz. 


-t 


-fan    fj"2A_e^_hi.  )  _  -fan"  /  m3A  e+  -  h+ 


Z   ft: 


Z/L* 


(HB) 


+an    I    W±Le+LzJ*\  _ 


-/-    -ran 


z  /^ 


■/an'  /^  e<  -  h* 


where 

d 


•4-        — 


=  /(?.- 

?,)^ 

Cu- 

?.r 

=       V(?3- 

-  ^f- 

>-cv 

-i,r 

=  ^T 

-*,r 

■+■  U. 

*--?3)r 

41 


y^,  -  ^fi-  o.-^)x 


(12B) 
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in    which 


ml3_  = 


h-    7, 


/r? 


*3 


m 


3* 


h  -    7* 


>"«,  = 


7v-  ?* 


7,  -  ? 


±_ 


(13B) 


and 


n 


=  V6r-Sfer  *-    (y-  ?k)^-t-    z"-'     ,  k=  i,z,3,4- 


(14B) 


=         Z*    +-     £x-*>k) 


fe=  /,<Ly3,4- 


(15B) 


h^  =  (y-7j(*-^ 


/e  =  /,  ^  3,  *- 


(16B) 


In  actually  evaluating  these  expressions  d>     and  <^>  cause 
no  trouble.   They  become  infinite  on  the  edges  of  the  quadrilateral, 
but  in  practice  they  are  never  evaluated  there.   The  component 
q>       requires  special  handling  in  certain  cases.   As   ^  -*  o  , 
<±)  -p.   q   if  the  point  P  is  approaching  a  point  in  the  plane  out- 
side the  boundaries  of  the  quadrilateral.  If  P  approaches  a 


point  within  the  quadrilateral  cj>  -*■   <27T(sgn  z)  as 


H  -#■  o 


These  facts  may  be  verified  from  equation  (11B).   In  the  course 
of  this  method  of  flow  calculation  it  is  required  to  evaluate 
<fi     at  points  in  the  plane  of  the  quadrilateral  elements.   In 
particular,  the  centroid  of  each  element  is  in  the  plane  of  that 
element  and  within  the  quadrilateral.   At  such  a  point  <p       should 
equal  Z1T    ,  since  the  case  of  interest  is  that  for  which  z:  — *  <=> 
through  positive  values,  rather  than  through  negative  values, 
i.e. ,  the  field  point  approaches  the  caisson  surface  from  the  ex- 
terior flow  field  rather  than  from  the  interior  of  the  caisson. 


142 


It  may  also  be  required  to  evaluate  induced  velocity  components 
at  points  in  the  plane  of  a  quadrilateral  outside  the  bound- 
aries of  the  quadrilateral,  for  example,  at  the  null  points 
(centroids)  of  other  elements  if  the  body  surface  has  a  flat 
region.   Points  in  the  plane  of  a  quadrilateral  element  should 
have  z  =  0,  but,  because  of  round-off  error,  they  may  have  have 
small  values  of  z  with  either  sign.   Thus,  for  points  inside 
the  quadrilateral,  equation  (11B)  may  give  ZTT    with  either  sign. 
To  avoid  this  error,  the  absolute  value  of  z  is  tested  before 
velocities  are  computed,  and  if  it  is  less  than  some  small  pre- 
scribed number,  which  is  nevertheless  large  compared  to  the 
expected  round-off  error,  it  is  set  equal  to  plus  zero  and  each 
inverse  tangent  of  equation  (11B)  is  set  equal  to   7T/2  with 
the  sign  of  the  numerator  of  its  argument.   This  gives  cp>  —  O 
for  points  outside  the  quadrilateral  and  cft>   =-  ztt       for  points 
inside  the  quadrilateral.   Another  situation  that  may  cause 
trouble  in  the  computing  machine  is  when  the  slope  of  a  side 
of  the  quadrilateral  is  infinite,  i.e. ,  when  a  side  is  parallel 
to  the  y-axis.   It  is  evident  from  equation  (11B)  that  in  this 
case  the  two  inverse  tangents  corresponding  to  that  side  cancel 
each  other.   To  avoid  difficulties  each  of  the  quantities 
(€*-"£,),  (  T3  -  7t   )  ,  (7+-%),    and  (  J  f    -  rf  )  are  tested 
to  determine  whether  they  are  zero,  and  if  any  one  of  them  is 
zero,  the  two  inverse  tangents  corresponding  to  that  side  are 
set  equal  to  zero.   Finally,  it  should  be  mentioned  that  the 
inverse  tangents  in  equation  (11B)  are  evaluated  in  the  normal 
range  -7T/2  to  -t-TT/2.       It  is  tempting  to  combine  some  of  the 
inverse  tangents  in  this  equation  using  the  tangent  addition 
law,  but  if  this  is  done,  great  care  must  be  exercised  with 


143 


regard  to  the  range  in  which  the  resulting  inverse  tangents 
should  be  evaluated. 

When  the  distance  R  is  large  it  is  possible  to  approximate 
integrals  in  Eq.(lB)  and  (6B-8B)  by  simply  multiplying  the  integrand 
by  the  area  AS.  Thus,  in  order  to  get  some  indication  of  how  large 
R  must  be  in  order  for  this  approximation  to  be  valid  Figure   (IB). 
The  results  shown  here  indicate  that  in  the  case  of  either  the 
induced  velocity  or  the  potential  contributed  by  the  panels  of  area 
AS  it  is  necessary  that  R2  2^A  S   in  order  for  the  approximate 
expression  to  become  valid. 
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APPENDIX  C:   ALTERNATE  FORM  OF  THE  GREEN'S  FUNCTION  CONVENIENT  FOR 
NUMERICAL  EVALUATION 

Except  when  the  point  x,y,z  lies       close  to  the  source 
located  at  £j",^,J  the  series  form  of  the  Green's  function  given  by 
Eq.(3.4  )  is  used  for  computer  numerical  evaluation.   The  series 
appears  to  converge  fairly  rapidly  and,  therefore,  the  numerical 
evaluation  is  quite  efficient.   However,  when  the  point  x,y,z  lies 
close  to  Jf,^,^   the  series  form  cannot  be  used  because  of  the 
singular  nature  of  the  Bessel  functions  yo(<2A) and  fC(&*-)as  J~L~Jf*   °  • 
In  such  cases,  i.e.,  when  (  Q/l)  becomes  less  than  some  arbitrarily 
small  number,  the  form  of  the  Green's  function  specified  by  Eq.(3.3  ) 
must  be  used. 

However,  Eq.(3.3)  may  not  be  evaluated  numerically  in  a  direct, 
straightforward  manner  because  the  1/R  term  is  singular  as  /^-*-  O 
and  the  integrand  of  the  infinite  integral  is  singular  at  jj.  =  a. 
The  evaluation  of  the  1/R  term  through  numerical  integration  was 
discussed  in  Appendix  B.   In  this  appendix  an  efficient  method  for 
the  numerical  evaluation  of  the  infinite  integral  is  discussed. 

The  integral  form  of  the  Green's  function  is  given  by 


(i-c) 
+  z  pv  Ccv+yie"* cos/?M7+h)i  coshfrfy+h)!  Jl (**)  J^ 

+  l  i?77Va*-vz)  cnshLaa  +  fol  c.oshfacUi-h)]  Jl(a*) 
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^  '  -^  S/nh  uj/i    -  >-  cosh  Mh 


where     /R  =  ^  (X-^  +   (y-?>*  +  &-~S)* 

V  =  CTzg  _   q  +anh(ah) 
9 

The  integrand  of  the  infinite  integral  can  be  simplified  by  use  of 

the  relationships,  • 

sjr?h(x)    ~  cosh  C*)     ^    ^  <S 

for  x  larger  than  about  4.0.   Thus,  the  interval  may  be  broken  into 
two  parts.   Using,  for  brevity 

■h/^ih 
Eq.(l-C)  may  be  written  as 

2  Rv.  \  (     )  da   =  Z  RV.    I    /         )  dyu 

r°°  .  s  (3-C) 

+  RV.f  (A4±2\eM<**»    TJM^    d* 

where  u     >  4,o/h}    .   The  term  (m+V)  /(M-y )  may  be  written 
so  that  the  infinite  integral  becomes 

z  p.v  I  (      )  djx.   =    2  pk  J  (       )   q^ 

f  /-'k'  [  _^>L_  e*"6"^  J2CMK)  da 


(4-C) 


^,    X" 


This  may  be  further  rearranged  by  use  of  the  relationship  given 
for  example,  by  Gradsheteyn  and  Ryzhik  (21). 
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,<*> 


J_   =.   f  e^C^V       ZTo(MA-)      <JM  (5-C) 

R"  I 

where  R"  represents  an  image  source  with  respect  to  the  free  sur- 
face 


y  '  (6-c) 


The  infinite  integral  may  now  be  written, 


o 

,.,.,,  (7-C) 


Z  RV.[(       )  4m    =  -L-„       +      2.  rvJJ         )    c/ 

The  remaining  infinite  integral  can  now  be  expressed  in  terms  of 
a  finite  integral  by  use  of  the  relationship  (See  Abramowitz  and 
Stegun(22) ) . 

Jo(MrL)    =  M  I     e  C?fe  (8-C) 

Using  Eq.(8-C),  the  last  term  may,  therefore,  be  written  as 

°°   TT 


Now,  let  t  =  >i-  V   and  dt  =  dj^  in  Eq.(9-C): 

.1  ^y-^        J.  f"0  4*  •  &i>      -£ ^-<2_'      oM(  10.C) 


or,  simplifying  further  by  use  of  the  substitution  S  =  t/Cu^-P) 


", 


5  =/  o 
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The  exponential  integral  as  defined  by  Abramowitz  and  Stegun  is 
so,  the  above  result  may  be  written  as 

OO 

Rvf-gi-r  eM(**V    Jit**.)    d^ 

\  .  (13-C) 

IT  J 


With  this,  the  complete  expression  for  G  may  be  written 

&  K'         F?"  '  J0  Cosh/lh     (y«UnhM     -y)  (14-C) 

<?       '  j£6«/o  d^  +  £?  e  E.L-Cyi-i+iA.coseXy,-*! da 

a  o 

This  form  of  the  Green's  function  involves  no  infinite  integrals 
and,  therefore,  can  be  evaluated  numerically  fairly  efficiently. 

The  integrand  of  the  first  integral  is,  however,  singular 
when  yU.  =  a  so  a  special  procedure  must  be  followed  in  evaluation 
of  that  integral .   Let 

ft*)  =  z(M+y)eM/l ca5h[M(i+h)~l  cosh[M(y-hh)l  J~0(^^)       (15_c) 

Cosh(<uh) 


The  integral  in  question  can  be  written 
U,  r2LQ 


za 


r   '  r  f 

Pyf    f<")f      cJm       _.    /    -F(m)    -    i(a)     X  +  f(a)\ 


hyuh-y 


a 


-h 


<"  -hxnh/xh  -  >> 
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The  first  and  third  integrals  in  Eq.(16-C)  are  not  singular  and 
will  cause  no  difficulty  in  numerical  integration  provided  the 
integrand  is  not  evaluated  exactly  at  the  midpoint  where  yu  =  a. 
The  integration  procedure  used,  therefore,  is  the  "3/8  rule"  which 
evaluates  the  integrand  at  end  points  and  two  intermediate  points. 
Using  an  odd  number  of  subdivisions  then  the  integrand  is  never 
evaluated  at  u  =  a. 

The  second  integral  in  Eq.(16-C)  is  singular  and  a  special 


procedure  must  be  followed.   The  integral  may  first  be  written  as 

(17-C) 


where 


M   ianh  Mh  -  V 


Uo  =    ah 


u 


U  -tanh  u~  u0  -hxnh  ut 


which  indicates  clearly  the  singular  nature  of  the  integrand  at 
IX-    UL0  .       The  interval  may  be  broken  into  three  parts  as 


f 


zu< 


du 


u0-e 


Ju 


J    U  far?/?  U  -  U0  tanh  U0      -  j   a   tank  U-  U0  -hanh  u\ 


(18-C) 


o 

rzua 
"*"  I u  tanh  u'-u0  fanh ul     *"       U  tanh  u       u0  +anhu0 

and  the  first  and  last  integral  integrated  in  a  straightforward 
manner  by  numerical  integration.   The  second  integral  is  evaluated 
by  first  expanding  the  denominator  in  a  Taylor  series  about 
and  then  dividing  1.0  by  this  series.   The  result  of  this  is 

{?**"*  use. 

JU  tanh  U-U0  -runhUo      '  \  I  Qt(U-  U0)  Q^  +   ~Q,(~Qf  ~a;j(U~U°)"']du 


UD-6 


u   -G 
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a,  a         a,1 


+ 


>-(QL  -    «iU 


a,  \  a," 


)A....J 


dx 


-e 


(19-C) 


where 


a 


a2  = 


hknh(ah)  +    Q-h   sechL(ah) 
sech^ah)  ( I-  ah   +cmh(ah) 


Thus,  the  complete  integral,  Eq.(17-C)  may  be  written 


.za 


du 


Uo~£ 


m    -ha n huh   -  y 


All 


a    tank  U~  U*   ianh  Uc 


zu0 


r 


■+■ 


du 


U   fanh  u-U0  fanh  Ur 


-  2<£ 


U0+£ 


(20-C) 

sechz(ah)  ( I  -  ah   ignh  ah) 
(ianh  (ah)  +  ah    scch^ah))2- 


where  £      denotes  some  small  value  between  0  and  QH .   In  the  numer- 
ical evaluation,  £  has  been  set  to  0.10. 

The  final  form  of  the  Green's  function  in  the  form  appropriate 
for  numerical  evaluation  may  now  be  written  as 

y 


G-M 


I 

R'    +     R 


za  u, 

R"      J  ^  fanh  juh  -  y     ^      J  M  fanh  Mh 


+ 


ah-e 
f(a)f 


o«  -f-anh^u  -ah  fanh  (ah) 


zah 


dM 


M  -f-anhsU-  ah-fanhmh) 


ah+k 


-  ze 


5e 


cti~{ah)    Cl-ali    fnnhfah)')     f(a) 


(21-C) 


(fanh  (ah)  +   ah  SecJfcah)  ) 
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+  3jLf<E?>iWLeose>     Sl-Cy+l  +  tAcoseXurV)]    do 

•    27r(a^-  \>*-)  coshLaCM-h)]   cosh[Q-C^-»-h)J    3~0ca/L) 
alh   -  v*h    +     v 

where 

r    .  _  z(^*y)eM    coshL*a+h)i  cosht^ty+h)]    Jiu^) 
t(M)  -  cosA^ 

The  derivatives  of  G(x,y,z; f  ,^  , J )  follow  from  Eq.(21-C)as 


o  ^ 


illi        «=  r- 


SgcH_oJi   C I-  ah  -fan/?  ah  )__    /Va;     (*-J.) 

(22-C) 


Q-f-anh  ah    +  ah    Sechfah  ) 


-f- 


(X-5)/^  <2^^?;       iSJ^S2   ^ 


rr      /t    J.  ^ 


(y+}+  I  re  cose  j  J 


I   z-rrcLLQSvll^  coskL'trM  «*hb<V+M  ^^  ^"^ 
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where 


Zajl  (M-hy)e^ coshfuO-t-h)]  COshL^cy-hh_)l     J^Cmsl) 


TW~-  cash(Mh)  "- 


za  .  u, 

+  I^-A~    f(a)       &      +f         J'f*  LC/M~- 


za 


i      J  m.  fa.nh.yu-  aft  -tank  ah  v      J  suJanh/i-  dh  -fanh(ah) 


Qh-t£ 


_  ze    sGah^ah  ( ' / '-  ah  iunh  ah  )     f^a)  (23-c) 

(■fdnh  ah   -h  ah  se.ali'ah  ) 


7r 
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ty-f-ii-L/z  cos  e)  U.    _ 

e      :  1  de 


+  [    *?m  (Q*--y*-)  cosh  Lac  ^3)]  smhrachi-y)l      T/n     , 


a'-h    -  -y^h    -h    >> 


where 


-f(M)  rr    ZM  (M^  y) )  e  Mh  CoshLP-Ci+foJSinhLuLH+h)]  Jra  (M/t) 
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and    finally, 


ak-e  zah 

-(z-1)  fU)\  U(  +anh£_  ahfanhah  "  &~Vf®JjT757^-ah-h3'ih  ah 


o  -■"* 


+  ze  s<s<zhzah    f  I  -  ah  +anh  ah)    f'(a)   (Z-y) 
(■fan/?  ah  -t  ah  sech  &h~) 


(24-c) 
Si. 


+■  i  *y(z^)&J>>e"»^  +  l'ca's6>E,£-e!/*?'tAcOIe>f<<.->'>] 


_    _e — L — -  /  do 
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